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ABSTRACT 


Air-sea  exchanges  of  heat  and  momentum  are  important  elements  in  understanding  and 
skillfully  predicting  tropical  cyclone  intensity,  but  the  magnitude  of  the  corresponding 
wind-speed  dependent  bulk  exchange  coefficients  is  largely  unknown  at  major  hurricane 
wind  speeds  greater  than  50  m  s'1.  Since  direct  turbulent  flux  measurements  in  these 
conditions  are  extremely  difficult,  the  momentum  and  enthalpy  fluxes  were  alternatively 
deduced  via  axisymmetric  angular  momentum  and  total  energy  budgets.  A 
comprehensive  error  analysis  was  performed  using  both  idealized  numerical  simulations 
to  quantify  and  mitigate  potentially  significant  uncertainties  resulting  from  unresolved 
budget  terms  and  observational  errors.  An  analysis  of  six  missions  from  the  2003 
CBLAST  field  program  in  major  hurricanes  Fabian  and  Isabel  was  conducted  using  a 
new  variational  technique.  This  analysis  indicates  a  near-surface  mean  drag  coefficient 
(Cd)  of  2.4xl0"3  with  a  46%  standard  deviation  and  a  mean  enthalpy  coefficient  (Ck)  of 
1.0x10"  with  a  40%  standard  deviation  for  wind  speeds  between  52  and  72  m  s'  .  These 
are  the  first  known  estimates  of  Ck  and  the  ratio  of  enthalpy  to  drag  coefficient  (Cr/Cd) 
in  major  hurricanes.  The  results  suggest  that  there  is  no  significant  change  in  the 
magnitude  of  the  bulk  exchange  coefficients  estimated  at  minimal  hurricane  wind  speeds, 
and  the  ratio  Cr/Cd  is  likely  less  than  0.75  for  wind  speeds  greater  than  50  m  s'1. 
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I. 


INTRODUCTION 


Heat,  moisture,  and  momentum  exchange  at  the  air-sea  interface  are  primary 
processes  in  hurricane  intensification  and  maintenance  with  important  theoretical  and 
practical  implications  (Kleinschmidt  1951;  Malkus  and  Riehl  1960;  Ooyama  1969; 
Emanuel  1986,  1995b).  Wind-speed  dependent  bulk  aerodynamic  formulas  often  are  used 
to  represent  the  turbulent  fluxes  associated  with  air-sea  interactions,  relying  on  drag  (Cd), 
moisture  (Ce),  and  enthalpy  (CK)  exchange  coefficients.  However,  as  Ooyama  (1969) 
stated,  “Unfortunately,  there  is  little  information  on  Ce  under  hurricane  conditions,  other 
than  the  semispeculative  guess  that  the  exchange  coefficients  of  latent  heat,  sensible  heat 
and  momentum  are  probably  of  the  same  magnitude.”  Nearly  40  years  after  that 
pioneering  study,  enthalpy  and  momentum  exchange  coefficients  are  still  largely 
unknown  at  major  hurricane  wind  speeds  (>50  m  s'1,  equivalent  to  category  three  and 
higher  on  the  Saffir-Simpson  scale).  There  are  known  sensitivities  of  theoretical  and 
numerical  models  of  major  hurricanes  to  the  magnitude  of  these  exchange  coefficients.  It 
is  therefore  important  to  improve  our  understanding  and  forecasts  of  these  phenomena, 
given  the  high  impact  that  their  damaging  winds  and  storm  surge  can  have  on  coastal 
populations,  global  economics,  marine,  and  Naval  operations. 

It  is  very  difficult  to  directly  measure  hurricane  surface  layer  fluxes  over  the 
ocean  due  to  the  extreme  conditions  and  challenging  deployment  of  instrumentation.  In 
situ  measurements  by  ships  or  manned  aircraft  are  very  hazardous,  and  fixed  sensors  are 
unlikely  to  be  in  the  correct  location  or  robust  enough  to  obtain  the  needed 
measurements.  Due  to  these  hazards  and  difficulties  in  obtaining  direct  air-sea  flux 
measurements,  only  fluxes  at  higher  altitudes  and  minimal  hurricane  force  winds  have 
been  obtained  (Black  et  al.  2007;  French  et  al.  2007;  Drennan  et  al.  2007;  Zhang  et  al. 
2008).  Measurements  of  the  drag  and  enthalpy  exchange  coefficients  at  high  wind  speeds 
also  have  been  made  in  laboratory  experiments  (Donelan  2004;  Haus  et  al.  2010),  but  are 
limited  to  wind  speeds  less  than  -50  m  s'1.  Global  positioning  system  (GPS) 
dropwindsonde  profiles  were  used  to  calculate  flux  profiles  and  deduce  the  drag 
coefficient  at  the  highest  wind  speeds  yet  (Powell  et  al.  2003;  Vickery  et  al.  2009),  but 
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this  technique  has  not  been  attempted  for  retrieving  the  enthalpy  exchange  coefficient. 
Indirect  retrievals  of  the  drag  coefficient  in  hurricanes  have  been  conducted  using  ocean 
measurements  (Shay  and  Jacob  2006;  Jarosz  et  al.  2007),  but  this  methodology  is  not 
applicable  for  deducing  enthalpy  exchange  coefficients. 

An  indirect  approach  to  retrieving  momentum  exchange  using  an  angular 
momentum  budget  was  originally  proposed  over  fifty  years  ago  (Palmen  and  Riehl  1957). 
In  this  formulation,  the  surface  stress  is  obtained  by  measured  transports  of  absolute 
angular  momentum  in  an  axisymmetric  cylindrical  coordinate  system.  If  the  flux  of  this 
quantity  is  known  at  the  top  and  sides  of  a  prescribed  control  volume,  the  flux  at  the  air- 
sea  interface  can  be  obtained  via  residual.  A  similar  budget  can  be  derived  for  total 
energy  to  deduce  the  enthalpy  exchange  coefficient.  The  central  focus  of  this  research 
was  the  application  of  this  budget  method  to  major  hurricanes  using  recent  observations. 
After  Palmen  and  Riehl  (1957)  introduced  the  budget  methodology,  additional  studies 
used  the  angular  momentum  budget  to  deduce  the  drag  coefficient  with  improved  datasets 
(Miller  1962,  1964;  Hawkins  and  Rubsam  1968;  Hawkins  and  Imbembo  1976).  The 
quality  and  density  of  hurricane  observations  has  improved  significantly  since  that  time, 
as  well  as  the  analysis  techniques  used  to  composite  the  data. 

In  practical  application  for  this  dissertation  research,  some  of  the  budget  terms  were 
very  difficult  to  calculate  with  observational  data,  and  known  uncertainties  also  must  be 
acknowledged.  To  gain  insight  into  the  relative  importance  of  the  various  sources  of  error  in 
the  energy/momentum  budget  method,  the  sensitivity  to  errors  in  unresolved  budget  terms, 
sea-surface  temperature  (SST),  center  placement,  gridding  method,  and  size  of  control 
volume  were  examined  systematically  using  numerical  model  data  where  the  surface  fluxes 
were  known.  These  results  were  then  utilized  to  determine  a  quantitative  estimate  of  the 
confidence  level  of  the  magnitude  of  the  retrieved  bulk  exchange  coefficients.  A  high- 
resolution  dataset  collected  in  Hurricanes  Fabian  and  Isabel  (2003)  as  part  of  the  Coupled 
Boundary  Layers  Air-Sea  Transfer  (CBLAST)  experiment  (Black  et  al.  2007)  was  used  to 
apply  this  energy/momentum  budget  method  to  real  tropical  cyclones.  Both  tropical  cyclones 
(TCs)  achieved  category  four  intensity  and  the  data  collected  present  a  significant  advance  in 
the  measurement  of  major  hurricanes.  Although  comprehensive  analyses  of  the  structure  of 
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Hurricane  Isabel  have  been  performed  elsewhere  using  dropsonde,  in  situ  flight  level, 
Doppler  radar,  and  satellite  data  (Montgomery  et  al.  2006;  Aberson  et  al.  2006;  Bell  2006; 
Bell  and  Montgomery  2008),  one  unique  aspect  of  the  dataset  has  yet  to  be  fully  utilized.  A 
series  of  rapid  dropsonde  releases  in  the  eyewall  region  of  these  two  hurricanes  were 
performed  to  accurately  construct  control  volume  composites  for  use  in  the  budget  method. 
These  dropsonde  “sequences”  are  used  in  conjunction  with  Doppler  radar,  in  situ,  and 
radiometer  data  to  construct  composite  analyses  for  the  air-sea  exchange  calculations 
presented  in  this  dissertation. 

To  deduce  momentum  and  energy  fluxes  from  these  sequences,  the  data  must  be 
gridded  in  an  axisymmetric  coordinate  system  moving  with  the  tropical  cyclone.  A  simple 
objective  analysis  scheme  (Barnes  1973)  was  used  in  the  previous  Isabel  studies  to  obtain 
the  kinematic  and  thermodynamic  structure.  Recent  advances  in  data  assimilation 
techniques  will  be  used  here  to  improve  this  analysis  by  deriving  a  variational  procedure 
that  provides  a  maximum  likelihood  estimate  of  the  gridded  structure  given  estimates  of 
background  and  observational  errors.  Additionally,  this  system  allows  for  the  incorporation 
of  Doppler  radar  data,  which  increases  the  spatial  coverage  and  improves  the  kinematic 
measurements  in  the  analysis.  This  integrated  variational  technique  also  can  provide  a 
unique  framework  for  future  theoretical  TC  studies  and  numerical  model  initialization. 

To  provide  a  broad  scientific  context  for  the  current  research,  Chapter  II  is  a  review 
of  relevant  previous  studies  on  air-sea  fluxes  in  tropical  cyclones.  A  derivation  of  the 
equations  for  total  energy  and  angular  momentum  conservation  in  an  axisymmetric, 
cylindrical  coordinate  system  necessary  for  calculations  using  the  budget  method  is  provided 
in  Chapter  III,  and  this  is  followed  by  a  comprehensive  error  analysis.  This  error  analysis 
includes  results  from  numerical  simulations  of  a  hurricane  vortex,  wherein  the  sensitivities  to 
simulated  observational  deficiencies  are  explored.  The  variational  methodology  for  deducing 
the  axisymmetric  tropical  cyclone  stmcture  from  real  observations  is  described  in  Chapter 
IV.  The  results  of  the  analysis  using  the  described  methods  are  presented  in  Chapter  V.  The 
concluding  chapter  summarizes  the  key  findings  of  this  study  and  discusses  potential 
implications  for  both  theoretical  understanding  of  hurricane  intensity  and  numerical  weather. 
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II.  TROPICAL  CYCLONE  SURFACE  FLUXES 


A.  THEORETICAL  AND  PRACTICAL  IMPLICATIONS 


Reynolds  decomposition  of  the  momentum  and  thermodynamic  equations  leads  to 
the  well-known  turbulence  closure  problem,  in  which  the  number  of  unknown  variables 
exceeds  the  number  of  equations  for  the  prediction  of  atmospheric  flow  (Stull  1988).  In 
the  surface  layer,  where  the  turbulent  fluxes  are  nearly  constant,  the  fluxes  are  often 
parameterized  by  bulk  aerodynamic  formulae  that  require  wind- speed  dependent 
exchange  coefficients  for  momentum,  heat,  and  moisture.  In  a  tropical  cyclone,  the 
surface  stress  can  be  represented  in  a  cylindrical  coordinate  system  in  which  the 
tangential  (swirling)  wind  is  the  dominant  component  as 

Ue  =  ~pwV  =  CDp  I  uh  I  v  (?  n 


where  v  denotes  the  tangential  wind,  w  the  vertical  wind,  p  the  density,  I  uh  I  the 

horizontal  wind  speed,  Cd  the  drag  coefficient,  primes  denote  perturbation  quantities,  and 
the  overbar  is  a  Reynolds  averaging  operator.  In  the  axisymmetric  coordinate  system 
used  in  this  study,  the  Reynolds  averaging  operator  encompasses  both  azimuthal  and 
temporal  averaging.  A  similar  expression  for  the  enthalpy  flux  can  be  derived: 

F:k  =  F:T  +  F:„  =  PCP  W'T'  +  PFw'q'  =  C \p  I  U h  I  (**-*)  ^  ^ 

where  T  denotes  the  temperature,  cp  the  specific  heat  of  air,  L  the  latent  heat  of 
vaporization,  q  the  water  vapor,  k  the  enthalpy,  k*  the  saturation  enthalpy  at  the  sea 
surface,  and  Ck  the  bulk  enthalpy  exchange  coefficient  (Emanuel  1995).  Since  the  ocean 
surface  characteristics  change  significantly  with  increasing  wind  speed  (Black  et  al. 
1986),  a  simple  extrapolation  of  the  bulk  exchange  coefficient  magnitudes  derived  at  low 
wind  speeds  to  50  m  s'1  and  beyond  is  not  necessarily  justified.  Emanuel  (2003) 
presented  a  similarity  hypothesis  for  enthalpy  and  drag  coefficients  at  very  high  wind 
speeds  using  dimensional  analysis,  and  proposed  that  their  magnitudes  become 
independent  of  wind  speed  at  approximately  hurricane  force  and  higher  (>33  m  s'1).  That 
analysis  suggested  also  that  the  enthalpy  coefficient  should  vary  as  a  function  of 
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temperature,  such  that  the  ratio  Ck/Cd  would  be  approximately  unity  at  typical  near¬ 
surface  air  temperatures  in  major  hurricanes.  The  earlier  work  of  Ooyama  (1969), 
Rosenthal  (1971),  and  Emanuel  (1986,  1995)  suggests  a  strong  dependence  of  Ck/Cd  on 
intensity.  Therefore,  validating  or  rejecting  these  hypotheses  about  the  behavior  of  the 
bulk  exchange  coefficients  at  major  hurricane  wind  speeds  has  important  implications  for 
tropical  cyclone  intensity  theory  and  numerical  weather  prediction  (NWP). 

The  importance  of  air-sea  energy  and  momentum  exchange  in  hurricanes  was  first 
introduced  by  Kleinschmidt  (1951)  in  the  form  of  an  energy  balance  between  the  amount 
of  sensible  and  latent  heat  a  TC  can  extract  from  the  ocean  surface  and  the  momentum 
lost  to  the  sea  through  frictional  dissipation.  Malkus  and  Riehl  (1960)  presented  a 
potential  intensity  (PI)  theory  based  on  the  concept  of  a  steady-state  TC  at  maximum 
intensity  that  maintains  this  energy  balance,  and  they  derived  an  expression  for  the 
maximum  tangential  wind  of  a  TC  given  the  a  priori  SST,  air  temperature,  and  bulk  heat 
and  momentum  exchange  coefficients.  While  they  acknowledged  the  importance  of 
latent  heat  fluxes  from  the  sea  surface,  the  bulk  moisture  exchange  coefficient  was  not 
explicitly  included  in  their  equation  for  maximum  tangential  wind.  Ooyama  (1969)  and 
Rosenthal  (1971)  were  the  first  to  explicitly  link  the  flux  of  water  vapor,  or  latent  heat, 
with  numerically  simulated  hurricane  intensity.  Emanuel  (1986)  expanded  these  concepts 
and  derived  an  expression  for  maximum  intensity  assuming  a  steady-state  balance 
between  the  amount  of  energy  obtained  from  the  ocean  and  that  dissipated  by  friction  in 
the  boundary  layer,  which  is  similar  in  general  terms  to  that  presented  by  Kleinschmidt 
(1951)  and  Malkus  and  Riehl  (1960).  The  primary  parameters  that  govern  Emanuel’s  PI 
(E-PI)  are  the  sea-surface  temperature,  the  outflow  temperature,  boundary  layer  relative 
humidity,  and  the  ratio  of  the  bulk  enthalpy  (Ck)  and  momentum  (Cd)  exchange 
coefficients  (Ck/Cd).  Due  to  limited  high-resolution  surface  flux  measurements  in  the 
extreme  conditions  of  a  hurricane  eyewall,  exchange  coefficients  are  the  most  uncertain 
of  these  parameters  at  wind  speeds  above  category  one  (Black  et  al.  2007). 

Bell  (2006),  Montgomery  et  al.  (2006),  and  Bell  and  Montgomery  (2008)  carried 
out  an  extensive  test  of  the  E-PI  predictions  using  observations  collected  during 
Hurricane  Isabel  as  part  of  the  2003  CBLAST  field  program.  Since  little  guidance  existed 
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for  the  magnitude  of  the  bulk  exchange  coefficients  at  major  hurricane  wind  speeds  for 
the  main  potential  intensity  estimate,  the  ratio  of  bulk  enthalpy  and  momentum  exchange 
coefficients  was  assumed  to  be  unity  (Ck/Cd  =1).  The  range  of  E-PI  estimates  from  Bell 
and  Montgomery  (2008)  is  -30  m  s'1  by  changing  the  assumptions  in  the  calculations 
(Fig.  1,  note  that  all  figures  appear  at  the  end  of  the  subsections).  The  upper-bound 
estimate  neglects  entirely  the  ocean  cooling,  includes  dissipative  heating,  and  assumes 
that  Ck/Cd  =  1.  The  lower-bound  estimate  assumes  Ck/Cd  =  0.5  and  that  the  ocean 
cooling  effect  is  compensated  by  the  dissipative  heating  effect.  Accurate  determination  of 
the  bulk  exchange  coefficients  would  clearly  help  narrow  the  uncertainty  of  these 
estimates. 

Recent  research  (Smith  et  al.  2008;  Bryan  and  Rotunno  2009)  has  suggested  that 
the  primary  discrepancy  between  observed  and  simulated  axisymmetric  intensities  and 
the  E-PI  theoretical  limit  and  is  due  to  the  implicit  assumption  of  gradient  wind  balance 
in  the  boundary  layer.  In  three-dimensional  simulations  by  Montgomery  et  al.  (2010),  the 
maximum  intensity  has  a  much  lower  sensitivity  to  the  surface  drag  coefficient  than  in 
the  axisymmetric  simulations.  The  role  of  radial  transport  of  moist  entropy  across  the 
eye-eyewall  interface  in  two  dimensions  (Persing  and  Montgomery  2003)  and  three 
dimensions  (Eastin  et  al.  2002;  Cram  et  al.  2003)  is  also  not  currently  included  in  PI 
theory.  While  it  is  clear  that  a  revised  PI  theory  will  need  to  be  developed  to  address 
these  legitimate  limitations,  this  does  not  diminish  the  importance  of  estimating  surface 
fluxes  at  high  wind  speeds  as  part  of  accurate  determination  of  hurricane  maximum 
intensity. 

The  importance  of  understanding  of  air- sea  exchange  at  high  wind  speeds  is  not 
limited  to  the  theoretical  domain.  While  tropical  cyclone  track  forecasts  have  been  found 
to  be  relatively  insensitive  to  the  parameterization  of  surface  fluxes  of  enthalpy  and 
momentum,  the  ability  to  skillfully  predict  tropical  cyclone  intensity  and  structure  has 
been  found  to  be  dependent  on  meaningful  and  accurate  parameterization  of  these 
processes  (Davis  et  al.  2008).  Ooyama  (1969)  was  the  first  to  demonstrate  the  importance 
of  parameterized  air-sea  fluxes  for  growth  and  maintenance  of  a  tropical  cyclone  using  a 
simplified  three-layer  model  of  the  moist  atmosphere.  Rosenthal  (1971)  and  Rotunno  and 
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Emanuel  (1987)  further  demonstrated  this  importance  in  hurricane  models,  and  Emanuel 
(1995)  demonstrated  sensitivity  of  the  simulated  intensity  to  the  prescribed  surface 
exchange  coefficients  in  a  simplified  model  framework.  Subsequent  studies  by  Braun  and 
Tao  (2000),  Bao  et  al.  (2002),  and  Davis  et  al.  (2008)  using  full  physics  models  have  also 
shown  the  sensitivity  of  hurricane  intensity  forecasts  to  modifications  in  the  surface 
roughness  length  parameterization  used  to  calculate  heat  and  momentum  fluxes. 
Hurricane  simulations  using  the  Navy  Coupled  Ocean/ Atmosphere  Mesoscale  Prediction 
System  are  also  sensitive  to  the  surface  flux  and  sea  spray  parameterizations  (J.  Doyle 
2009,  personal  communication).  Despite  improvements  in  the  quality  of  the  intensity 
guidance  available  to  forecasters  at  the  National  Hurricane  Center,  relatively  little 
improvement  has  been  achieved  in  the  forecast  of  tropical  cyclone  intensity  since  1990 
(Rappaport  et  al.  2009).  Reducing  errors  in  the  numerical  weather  model  representation 
of  air-sea  interaction  in  major  tropical  cyclones  would  likely  contribute  to  improved 
intensity  forecasts  issued  by  the  Joint  Typhoon  Warning  Center  and  the  National 
Hurricane  Center. 
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Figure  1. 
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Theoretically  predicted  azimuthal  mean  Vmax  at  the  boundary  layer  top  for 
varying  outflow  temperature  and  near-core  SST  with  a  constant  RH  =  80% 
(dashed  line).  ‘X’  indicates  the  primary  potential  intensity  estimate  for  the 
observed  near  environment  around  Isabel  (using  27.5  °C  SST  near  eyewall 
region  associated  with  Fabian’s  wake),  which  yields  a  56.6  m  s'1  mean  Vmax- 
The  dark  solid  curve  represents  the  average  storm-relative  tangential  wind 
speed  at  the  top  of  the  boundary  layer  derived  from  the  dropwindsonde 
measurements.  The  shading  represents  the  6ms'1  standard  deviation  of  this 
mean  value.  This  a  priori  E-PI  estimate  assumes  Ck/Cd  =1.  From  Bell  and 
Montgomery  (2008). 
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B.  PREVIOUS  SURFACE  FLUX  OBSERVATIONS 


Surface  fluxes  are  typically  measured  by  one  of  three  methods-eddy  correlation, 
inertial  dissipation,  or  flux-profile  relationships.  Each  of  these  methods  has  different 
advantages  and  disadvantages  depending  on  the  measurement  apparatus  and  field 
conditions,  but  all  three  present  difficulties  in  the  TCBL.  At  major  hurricane  wind  speeds, 
the  main  challenge  for  these  types  of  measurements  is  placement  and  robustness  of  the 
sensors.  Given  that  the  spatial  and  temporal  occurrence  of  these  wind  speeds  is  very 
small  compared  to  the  ocean  basin  size  and  seasonal  timescales  in  which  hurricanes 
occur,  the  probability  of  a  fixed  sensor  encountering  them  is  very  low.  Additionally, 
operational  ocean  sensors  (such  as  buoys)  typically  are  not  designed  to  withstand  the 
extreme  conditions  found  in  the  open  ocean  at  the  base  of  a  major  hurricane  eye  wall,  and 
instrumented  ships  prudently  avoid  these  conditions  as  well.  Therefore,  laboratory 
experiments  and  observations  from  research  aircraft  deployed  in  tropical  cyclones  are 
currently  the  most  viable  choices  for  obtaining  information  about  surface  fluxes  at  high 
wind  speeds.  However,  low-altitude  flying  in  the  turbulent  boundary  layer  at  the  base  of 
an  eyewall  is  also  very  hazardous,  which  requires  that  remote  sensing  or  expendables 
(i.e.,  dropwindsondes)  be  used  to  obtain  measurements  in  this  region.  This  generally 
limits  the  use  of  the  eddy  correlation  or  inertial  dissipation  methods  except  in  special 
circumstances  where  it  is  safe  to  fly  in  the  boundary  layer  and  the  aircraft  is  equipped 
with  high-resolution  sensors.  Similarly,  the  flux-profile  method  was  not  viable  until  the 
development  of  the  modern  NCAR  GPS  dropwindsonde  (Hock  and  Franklin  1999, 
hereafter  “dropsonde”).  It  is  perhaps  not  surprising,  then,  that  the  first  attempts  to 
determine  surface  fluxes  in  tropical  cyclones  were  performed  by  indirect  retrievals 
similar  to  that  proposed  in  this  study. 

Palmen  and  Riehl  (1957)  reported  the  first  calculation  of  the  drag  coefficient  in  a 
tropical  cyclone  using  a  tangential  wind  budget  using  two  composite  wind  fields  created 
by  Hughes  (1952)  at  low  levels  and  Jordan  (1952)  at  upper  levels.  Hughes  (1952) 
analyzed  84  flights  from  28  storms  in  which  reconnaissance  was  available  during  1945- 
1947,  while  Jordan  (1952)  used  rawinsonde  observations  collected  from  islands  and 
coastal  stations  near  tropical  cyclones  during  the  period  1945-1951.  These  two 
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composites  were  combined  to  create  an  axisymmetric  tropical  cyclone  structure  with  one 
degree  latitude  horizontal  resolution  and  variable  vertical  resolution  for  use  in  the  budget 
calculations.  The  original  resolution  of  the  low-level  data  was  0.5  degrees  with  the  lowest 
reported  winds  at  500  ft  (152  m),  and  no  upper-air  data  were  obtained  within  -200  km  of 
the  center.  To  derive  the  surface  stress,  the  tangential  momentum  equation  was  integrated 
from  the  surface  to  the  top  of  the  inflow  layer,  assuming  steady- state  and  cylindrical 
axisymmetry,  which  yielded: 


*:8 


(2.3) 


where  denotes  the  absolute  vorticity  /  +  d(rv)/rdr ,  and  pO  and  pH  are  the  pressures  at 

the  surface  and  top  of  the  inflow  layer,  respectively.  Vertical  motion  was  derived  from 
the  mass  continuity  equation  with  a  presumed  zero  boundary  condition  at  the  surface. 
This  calculation  yielded  a  Cd  of  2.1x10'  for  the  innermost  radial  ring  with  a  surface  layer 
wind  speed  of  -26  m  s’1,  but  it  apparent  that  this  was  a  coarse  estimate  given  the 
extensive  averaging  used  to  create  the  dataset.  Subsequent  budgets  of  angular  momentum 
and  energy  were  calculated  in  that  study,  but  no  computation  of  an  enthalpy  exchange 
coefficient  was  attempted. 


This  budget  method  was  repeated  in  several  subsequent  studies  with  improved 
datasets.  Miller  (1962)  used  an  aircraft  dataset  collected  in  Hurricane  Helene  (1958)  to 
extend  the  Palmen  and  Riehl  (1957)  results  by  using  a  single  storm  as  opposed  to  a 
composite,  and  at  higher  surface  wind  speeds.  This  study  yielded  drag  coefficients 
ranging  from  2.4x10  to  3.2x10'  for  winds  from  30  to  40  m  s'  .  Aircraft  data  from 
Hurricane  Donna  (1960)  were  used  to  calculate  drag  coefficients  of  3.6x10'  and  4.19 
xlO'  at  36  and  52  m  s'  ,  respectively  by  Miller  (1964),  as  shown  in  Figure  2.  These 
results  indicated  that  drag  coefficient  continued  to  increase  with  wind  speed  as  storms 
reached  major  hurricane  intensity.  Hawkins  and  Rubsam  (1968)  and  Hawkins  and 
Imbembo  (1976)  continued  these  analyses  with  Hurricane  Hilda  (1964)  and  Hurricane 
Inez  (1966),  and  extended  the  wind  speed  behavior  of  the  drag  coefficient  to  nearly  70 
m  s'  as  shown  in  Figure  3.  To  match  the  maximum  drag  coefficient  of  4.6x10'  at  67 

m  s'1  found  using  the  Inez  dataset,  the  least-squares  curve  fit  was  changed  from  linear  to 
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quadratic.  Hawkins  and  Rubsam  (1968)  also  included  a  more  comprehensive  error 
analysis  of  the  budget  method,  calculating  the  surface  stress  using  both  the  integrated 
tangential  momentum  equation  (so-called  “vorticity”  method)  and  the  angular  momentum 
form,  which  is  similar  to  the  current  study.  These  two  methods  are  analytically 
equivalent,  but  produced  slightly  different  numerical  results.  Also  particularly  relevant  to 
the  current  study,  Hawkins  and  Rubsam  also  compared  the  results  obtained  when 
changing  the  top  of  the  integration  volume  from  the  top  of  the  inflow  layer  to  the  top  of 
the  atmosphere  (100  hPa).  The  differences  between  the  two  integration  methods  and 
different  boundary  conditions  suggested  an  error  in  the  drag  coefficient  estimates  of 
about  ±30%.  An  energy  budget  was  also  calculated,  but  the  Ch  and  Ce  values  were 
assumed  to  be  equal  in  magnitude  to  the  derived  Cd,  as  opposed  to  being  retrieved  from 
the  budget  itself.  In  summary,  none  of  these  studies  using  budget  methods  indicated  that 
the  drag  would  level  off  or  decrease  at  major  hurricane  wind  speed,  and  none  presented 
estimates  of  a  heat  or  enthalpy  exchange  coefficient.  Frank  (1984)  also  used  an  angular 
momentum  budget  to  estimate  the  drag  coefficient  from  a  composite  analysis  of 
Hurricane  Frederic  (1979).  Frank  estimated  a  ±50%  error  in  Co  using  this  method  by 
calculating  a  standard  deviation  of  estimates  from  different  radial  bands. 

Moss  and  Rosenthal  (1975)  used  the  same  aircraft  dataset  from  Hurricane  Inez 
(1966)  to  calculate  the  drag  coefficient  using  the  Deardorff  boundary  layer 
parameterization  (Deardorff  1972),  and  their  results  compared  very  favorably  to  those 
obtained  by  Hawkins  and  Imbembo  (1976)  using  the  angular  momentum  budget.  They 
also  provided  an  estimate  of  the  bulk  heat  exchange  coefficient  at  48  m  s'1,  but  they 
admitted  that  it  clearly  had  a  high  bias.  Moss  and  Rosenthal  also  directly  compared  these 
drag  coefficients  with  the  simple  Deacon’s  formula  (Cd  =  1.1x10'  +  4x10'  IV1)  and  a 
formula  proposed  by  Miller  (1969,  Co  =  lxlO'3  +  7x10  5I VI),  and  concluded  that  these 
equations  gave  values  that  were  too  low  at  high  wind  speeds,  but  too  high  at  low  wind 
speeds  for  Inez,  but  were  reasonable  for  the  dataset  from  Hurricane  Daisy  (1958). 
Alhough  not  specifically  conducted  in  the  TCBL,  Large  and  Pond  (1981,  1982) 
performed  eddy  correlation  and  inertial  dissipation  estimates  of  open  ocean  heat  and 
momentum  exchange  in  wind  speeds  up  to  ~25  m  s'1,  and  their  results  have  been 
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validated  by  subsequent  studies  (Smith  et  al.  1992).  These  measurements  added  further 
support  to  the  hypothesis  that  the  drag  coefficient  continues  to  increase  with  wind  speed 
up  to  25  m  s'1.  In  contrast,  Large  and  Pond  (1982)  showed  no  significant  wind  speed 
dependence  on  the  moisture  flux  in  this  wind  speed  range,  and  this  result  has  also  been 
further  validated  (DeCosmo  et  al.  1996). 

Powell  et  al.  (2003)  used  a  flux-profile  relationship  to  infer  the  drag  coefficient 
from  dropsonde  profiles  released  in  hurricane  eyewalls.  The  wind  speed  is  obtained  via 
the  displacement  of  the  dropsonde’ s  GPS  position  as  it  descends  toward  the  ocean  surface 
after  release  from  a  research  aircraft.  Although  the  dropsonde  is  a  Lagrangian 
measurement  as  it  travels  several  kilometers  downwind  as  it  falls,  the  wind  profiles  can 
be  treated  in  an  Eulerian  framework  by  assuming  the  azimuthal  variation  in  the  wind 
speed  is  small  compared  to  the  vertical  variation.  Fitting  the  vertical  profile  to  a  least- 
squares  logarithmic  line  yields  the  surface  roughness  (zo)  as  the  intercept  and  the  friction 
velocity  (u*)  as  the  slope  from  the  flux-profile  relationship  in  a  neutral  surface  layer.  The 
implied  drag  coefficient  from  these  results  are  plotted  with  estimates  from  the  earlier 
budget  studies  in  Figure  4.  Powell  et  al.  provided  the  first  indications  of  a  decrease  in  the 
drag  coefficient  from  maximum  values  around  -2.5x10'  at  30-40  m  s"  .  Additional 
dropsonde  profiles  have  been  collected  in  major  hurricanes  since  the  2003  study, 
including  some  used  in  this  research,  which  have  extended  the  wind  speed  dependence  to 
near  60  m  s'1  (Figure  5)  with  a  near  constant  drag  coefficient  even  at  this  intensity 
(Vickery  et  al  2009). 

Additional  evidence  that  the  drag  coefficient  does  not  increase  above  -35  m  s"1 
was  reported  by  Donelan  et  al.  (2004)  using  laboratory  tank  measurements.  The  eddy 
correlation  method  was  used  for  wind  speeds  up  to  26  m  s'1,  and  a  momentum  budget 
retrieval  was  used  for  wind  speeds  from  20-50  m  s'1.  These  tank  measurements  showed  a 
“saturation”  of  the  drag  coefficient  around  hurricane  force  wind  speed  (33  m  s'1),  and 
suggested  a  limiting  aerodynamic  roughness  of  the  surface  waves  above  these  speeds. 
The  momentum  budget  approach  used  in  the  laboratory  has  also  been  utilized  in  field 
experiments  with  ocean  current  measurements.  Shay  and  Jacob  (2006)  used  airborne 
expendable  current  profilers  released  into  Hurricane  Gilbert  (1988)  to  estimate  the 
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downward  kinetic  energy  flux  and  retrieve  the  drag  coefficient  up  to  ~40  ms'1,  which 
leveled  off  above  28  m  s'1  to  a  value  near  3.5x10  3  (dashed  line  in  Figure  4).  A  similar 
approach  was  used  by  Jarosz  et  al.  (2007)  to  deduce  the  drag  coefficient  from  acoustic 
Doppler  current  profilers  moored  in  the  path  of  Hurricane  Ivan  (2006).  A  momentum 
budget  using  this  dataset  indicated  a  quadratic  relationship  with  wind  speed,  with  a  peak 
value  near  2.5x10  3  at  ~33  m  s'1  that  then  decreased  to  1.5xl0~3  at  45  m  s'1. 

The  uncertainty  in  the  magnitude  of  the  drag  and  enthalpy  exchange  coefficients 
in  the  TCBL  was  one  of  the  factors  leading  to  the  deployment  of  the  NOAA  P-3  aircraft 
as  part  of  the  CBLAST  field  campaign  (Black  et  al.  2007).  Due  to  the  safety  hazards 
mentioned  previously,  the  aircraft  were  only  flown  in  the  clear  air  boundary  layer 
between  rainbands,  but  fortunately  favorable  conditions  for  turbulence  observations  were 
achieved  in  two  major  hurricanes.  French  et  al.  (2007)  reported  the  measurements  of 
momentum  flux  from  48  flux  calculation  legs  flown  within  400  m  of  the  surface  in 
Hurricanes  Isabel  and  Fabian  (2003).  These  results  provided  the  first  open-ocean  eddy 
correlation  measurements  in  the  TCBL  at  hurricane-force  wind  speeds.  Estimates  of  the 
drag  coefficient  versus  wind  speed  from  this  study  (Figure  6)  indicated  no  discemable 
dependence  on  speed  in  the  range  measured.  Although  the  French  et  al.  results  slightly 
differ  from  those  reported  by  Powell  et  al.  (2003)  and  Donelan  et  al.  (2004),  general 
agreement  exists  that  the  extrapolated  Large  and  Pond  (1981)  formula  results  in  too  high 
drag  coefficients.  Latent  heat  (Drennan  et  al.  2007)  and  enthalpy  (Zhang  et  al.  2008) 
fluxes  derived  from  eddy  correlation  measurements  during  CBLAST  are  consistent  with 
a  lack  of  dependence  of  the  drag  coefficient  on  wind  speed  as  suggested  in  other  studies. 
Recent  laboratory  research  also  indicates  that  the  enthalpy  exchange  coefficient  is  nearly 
independent  of  wind  speeds  between  13  and  40  m  s'1  (Haus  et  al.  2010,  Jeong  et  al. 
2010).  Enthalpy  exchange  coefficients  as  a  function  of  wind  speed  from  the  HEXOS, 
CBLAST,  and  laboratory  experiments  taken  from  Haus  et  al.  (2010)  indicate  some  scatter 
but  no  significant  trend  with  increasing  wind  speed  above  -10  ms'1  (Ligure  7). 

In  summary,  most  recent  research  has  indicated  a  steady  or  slightly  decreasing 
drag  coefficient  with  increasing  wind  speed  beyond  30  m  s'1,  and  also  near-constant 
enthalpy  exchange  coefficient  in  hurricane-force  winds.  The  limited  estimates  of  the  drag 
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coefficient  above  50  m  s'1  are  inconclusive  as  to  the  magnitude,  with  earlier  budget 
studies  suggesting  a  substantially  higher  drag  than  the  flux-profile  method  of  Powell  et  al. 
(2003)  and  Vickery  et  al.  (2009).  Although  no  known  eddy  correlation  measurements 
exist  at  these  wind  speeds,  the  laboratory  tank  momentum  budget  in  Donelan  et  al.  (2004) 
is  in  closer  agreement  with  the  studies  of  Powell  et  al.  and  Vickery  et  al.  The  magnitude 
of  the  enthalpy  exchange  at  major  hurricane  wind  speeds  is  even  more  uncertain.  While 
there  currently  is  no  evidence  to  suggest  that  the  wind  speed  independence  should  not 
continue  above  35  m  s’1,  photographs  of  the  sea  surface  in  these  conditions  depict  a 
different  character  to  the  sea  state  than  in  minimal  hurricane  winds  (Black  et  al.  1986). 
The  sea  state  at  Beaufort  category  19  (>50  m  s’1)  shown  in  Figure  8b  is  described  as  “low 
clouds,  spray,  and  foam  merge  into  large,  white  areas  frequently  referred  to  as  ‘white 
sheets’  by  reconnaissance  crews”  that  cover  the  entire  surface  at  these  wind  speeds. 
White  water  only  covers  about  30-40%  of  the  photograph  at  -30  m  s’1  (Figure  8a),  with 
narrow,  parallel  streaks  being  the  defining  characteristic.  It  must  be  noted  that 
parameterizing  the  complexity  of  the  air-sea  interaction  at  these  wind  speeds  with  10  m 
bulk  exchange  coefficient  may  be  an  over-simplification,  but  given  the  established 
theoretical  and  numerical  reliance  on  these  quantities  it  is  still  believed  to  be  a  useful 
endeavor. 
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Figure  2. 


Drag  coefficient  Cd  as  a  function  of  wind  speed  over  water.  From  Miller 
(1964). 
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DRAG  COEFFICIENT 


Figure  3. 


DRAG  COEFFICIENTS  FOR  HURRICANE  INEZ  VS.  A  BACKGROUND 
OF  OTHER  HURRICANE  AND  LOWER  SPEED  DETERMINATIONS 


WIND  SPEED  (M/SEC) 

Drag  coefficients  for  Hurricane  Inez  (1966)  from  Hawkins  and  Imbembo 
(1976)  plotted  with  values  for  Hurricanes  Hilda  and  Helene  and  some  lower 
speed  determinations. 
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Drag  coefficient  as  a  function  of  10  m  wind  speed  from  dropsonde  flux- 
profile  relationship  (open  symbols)  from  Powell  et  al.  (2003)  and  previous 
studies  (closed  symbols,  lines). 
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Variation  of  the  drag  coefficient  with  mean  wind  speed  at  10  m  altitude  near 
the  radius  of  maximum  wind  (RMW).  From  Vickery  et  al.  (2009). 
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Estimated  drag  coefficient  as  a  function  of  Um  for  the  48  flux  runs  from 
French  et  al.  (2007)  showing  the  binned  values  (circles)  and  the  95% 
confidence  interval  from  this  study  and  extrapolation  of  results  from  Large 
and  Pond  (1981)  and  Smith  (1980),  dotted  and  dashed-dotted,  respectively. 
Also  shown  are  results  from  Donelan  et  al.  (2004;  diamonds)  and  Powell  et 
al.  (2003;  squares)  to  42  m  s'1. 
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Figure  7. 
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Wind  speed  dependence  of  Ck  from  Haus  et  al.  (2010).  ASIST  laboratory 
results  (•)  and  CBLAST  (A)  measurements  are  shown  with  HEXOS  results 
(x).  After  binning  observations  by  wind  speed,  the  mean  and  95%  confidence 
intervals  as  determined  from  a  t-distribution  of  the  combined  HEXOS  and 
CBLAST  field  data  are  shown  in  black. 
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Figure  8.  Typical  appearance  of  the  sea  surface  in  hurricane  conditions  of  Beaufort 

category  (a)  11  and  (b)  19.  The  photographs  were  taken  in  (a)  Hurricane 
Eloise  at  an  altitude  of  312  m  at  2246  UTC  22  September  1975.  The  aircraft 
was  in  the  right-front  quadrant  of  the  storm  about  140  km  from  the  center. 
The  flight-level  wind  was  32.8  m  s'1  and  the  mean  and  sustained  20  m  winds 
were  27  and  30  m  s’1,  respectively,  (b)  Hurricane  David  at  an  altitude  of  454 
m  at  1104  UTC  30  August  1979.  The  aircraft  was  in  the  right-rear  quadrant 
about  45  km  from  David’s  center.  The  flight-level  wind  was  63.0  m  s’1,  and 
the  mean  and  sustained  winds  were  50  and  57  m  s’1,  respectively.  Mean  and 
sustained  winds  are  ten  and  one-minute  averages,  respectively.  From  Black  et 
al.  (1986). 
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III.  BUDGET  METHODOLOGY 


A.  CONSERVATION  OF  ANGULAR  MOMENTUM 


To  determine  the  bulk  momentum  exchange  coefficient,  the  flux  form  of  the 
azimuthally-averaged  tangential  momentum  equation  in  cylindrical  coordinates  is: 


d{pv)  d(pruv)  d(pwv)  ,  puv 
- + - r - +  — r - + - +  p.fll  =  Fe 


dt 


rdr 


dz 


(3-1) 


where  u  is  the  radial  wind,  v  the  tangential  wind,  w  the  vertical  wind,  p  the  density, /the 
Coriolis  parameter,  and  Fe  molecular  friction/diffusion.  Multiplying  Equation  (3.1)  by  the 
radius  r,  and  assuming  that  Feis  small  above  the  viscous  sublayer  yields: 


d_ 

dt 


(pw)+ 


—  (p«rv)+  puv 
dr 


+  —  (pwrv)+  {pruf}=  0 
dz 


(3.2) 


By  the  chain  rule,  the  term  in  brackets  reduces  to: 


dipujrvp 

dr 


+  puv  ■ 


d(rpu{rv)) 

rdr 


(3.3) 


Considering  the  final  term  in  braces,  the  axisymmetric  mass  continuity  equation  is: 

dp  |  d(pru)  |  d{pw)  _ 

dt  rdr  dz  (3.4) 

Assuming  an  /-plane  approximation,  the  Coriolis  parameter  is  a  constant  and  therefore 
multiplying  the  continuity  equation  by  r//2  yields: 

r2d{pf)  |  rd(pruf)  |  r2d(pwf)  =  Q 
2  dt  2 dr  2 dz  (3  31 


By  the  chain  rule,  the  second  term  can  be  expressed  as: 


rd(pruf) 
2  dr 


diprufir1)) 
2  rdr 


-pruf 


(3.6) 
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Therefore: 


r  _  d(pfr2)  +  d(pru(fr2))  +  djpwfr2) 
2  dt  2  rdr  2  dz 


(3.7) 


Substituting  (3.3)  and  (3.7)  into  (3.2)  yields: 

|[p(rv  +  fr2/2)]  +  J ~[pru  (rv  +  jr2  / 2)]  +  ^-[pw (rv  +  jr2 /l)]  = 

d(pM)  d(pruM)  d(pwM) 
dt  rdr  dz 


where  the  absolute  angular  momentum  is  defined  as  M  =  rv  +  —  fr 


2 


Each  variable  is  then  separated  into  a  mean  and  perturbation  quantity,  such  that 
v  =  v  +  V  ?  where  the  mean  is  both  an  azimuthal  and  temporal  average.  Reynolds 
averaging  and  making  the  Boussinesq  approximation  (e.g.,  neglecting  p' )  yields: 


d(pM)  d(pruM)  d(pwM)  d(rpu'M')  d(pw'M') 
dt  rdr  dz  rdr  dz 


The  Reynolds  stresses  are  then  defined  asr rg  =  -pu'v'  and  Tz0  =  -pw'v' ,  where  the 
subscripts  represent  the  wind  directions  of  the  covariances  (r  is  the  radial  direction,  0  is 
the  tangential  direction,  and  z  is  the  vertical  direction).  Expanding  the  last  two  terms 

d(rpu'M')  3(r2  pu'v')  d(r2  xrg) 
rdr  rdr  rdr 

d(pw'M')  d(rpw'v')  r9r,e 

dz  dz  dz  (3.10) 


yields  the  axisymmetric  angular  momentum  equation  in  flux  form: 

d(pM)  |  d(pruM)  |  d(pwM)  d(r2 Trfj)  rdrze  _  Q 

dt  rdr  dz  rdr  dz 


(3.11) 
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Integrating  over  a  control  volume  from  z\  to  zi  and  r\  to  r2  yields  the  integrated  flux  form 
of  the  conservation  of  angular  momentum  in  axisymmetric  cylindrical  coordinates: 


n 

j  Zi  J  r 


rclrdz 


+ 


d(pM) 

zi  J'i  dr 

p  p  djpruM -r2zre) p  p2|~ d(pwM  -  rxz6) 
•hi  Jrt  Jr,  J: 


dr 


dz 


dz  rdr  ■ 


d(pM) 


P  P  ’Inrdrdz 

J;i  Jri  dt 

+ J  2  2 nr  [i puM  —  r  Ti  fj  ]  l(  dz  -  J  2 nr  [ puM  -  r  ti  B  ]l,  dz 

Z\  Z\ 

+  f  ~  [  pwM  -  rxzd  ]  L  2  nrdr  -  J  2  [  pwM  -  rxzd  ]  L  2  nrdr  =  0 


(3.12) 


Rearranging  (3.12)  and  dividing  by  2m ,  the  radially-integrated  turbulent 
momentum  flux  at  z\  is  then: 

r2dr  = 

-r2  J  [puM  -  rxre  ]  [2  dz  +  r,  J  -  rxr9  ]  I,,  dz 

-  J  *  [  pwM  -  rxz0  ]  L  rdr  +  |  ‘  [pvvAf  ]  I  rdr  -  J  "  J  ‘  — ^ — -rdrdz 


dt 


(3.13) 


From  (3.13),  the  integrated  surface  stress  for  an  arbitrary  control  volume  may  be 
calculated.  A  schematic  of  a  hypothetical  control  volume  is  shown  in  Figure  9,  which 
illustrates  the  approximate  location  of  zi,  Z2,  O,  and  r2  for  an  idealized  hurricane  flow. 
Note  that  zi  does  not  necessarily  have  to  be  at  10  m  altitude,  but  could  be  anywhere  in  the 
surface  layer  assuming  the  fluxes  are  nearly  constant  in  that  layer. 

Assuming  the  vertical  stress  at  z\  can  be  represented  by  a  bulk 

formula, 7:0  -  pw  v  -  CDp  I  uh  I  v  ,  an(j  negiecting  the  radial  variation  of  Cd  over  the 
control  volume  yields  the  drag  coefficient: 
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(3.14) 


cd  =(-r2^  [puM  ~rtre\ l2  dz  +  r,\  [puM -nr9]\ri  dz 

-J/  [ pwM  -  rrze ]  l  rdr  +  J;  [pwM }  I.  rdr  -  J*  J;  ^^-rdrdz ]  / 

J  [P 1 «/, 1  v]  1,  /'2^' 

w  rl 


Since  the  tendency  and  additional  flux  terms  cannot  be  calculated  using  the 
current  dataset,  the  working  form  for  use  with  these  observations  is 


CD  =  (-r2  J  [puM  ]  lr2  dz  + 1]  J  [puM  ]  I,,  dz 
-  J"2  [pwM  ]  I  rdr  +  J'"  [  pwM  ]  I  rdi\  /  ( J':  [p  I  uh  I  v]  I.  r2di\  +  R 


(3.15) 


where  the  unresolved  residual  is 


R  =  (,-J"[rr,.e]lr2  dz~r^[nrd]\  *+J'J  [rrz9]  [2  rdr-\\'  — 
J  :  [p  I  uh  I  v]  L  rdr 


^rdrdz  )/ 


(3.16) 


Alternatively,  one  can  obtain  the  average  surface  stress  by  dividing  (3.13)  by  the 

rr2  2 

integrated  square  of  the  radius  r  dr ,  or  divide  by  the  mass-weighted  term 

rr2  2 

pr'  dr  to  obtain  the  average  friction  velocity  u*  in  the  control  volume.  The  residual 

Jr, 

terms  were  also  unresolved  in  the  previous  studies  utilizing  a  budget  methodology.  An 
accurate  determination  of  the  surface  stress,  friction  velocity,  or  drag  coefficient  from  this 
formulation  requires  minimal  errors  in  the  specification  of  the  axisymmetric  mass  and 
wind  fields,  strict  adherence  to  axisymmetric  mass  continuity,  and  a  small  magnitude  of 
the  unresolved  residual  terms.  Each  of  these  error  sources  is  examined  in  Chapter  III.D. 
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Figure  9.  Schematic  illustrating  hypothetical  control  volume  (black  dashed  line)  used 

for  the  budget  methodology.  A  simplified  secondary  circulation  (gray 
streamlines)  and  region  of  maximum  wind  (vmax)  are  shown  to  indicate  the 
control  volume  encompasses  the  eyewall  region. 


B. 


CONSERVATION  OF  ENERGY 


To  determine  the  bulk  moist  enthalpy  exchange  coefficient,  the  first  law  of 
thermodynamics  in  material  form  with  enthalpy  as  the  state  variable  is: 


> D(yCpT)  =  ^  +  £v2r  +  qr  +  qf  +  ql 

Dt  Dt  L 


(3.17) 


where  cp  is  the  specific  heat  at  constant  pressure,  T  is  the  temperature,  k  is  the  thermal 
conductivity,  QR  is  the  diabatic  heat  exchange  due  to  radiative  transfer,  Qr  is  the  heat 
arising  from  frictional  dissipation  of  kinetic  energy,  and  Ql  is  the  latent  heat  release 
associated  with  phase  changes  of  water,  with  the  additional  terms  defined  above  in 
Chapter  III.A.  The  latent  heat  release  is  given  by: 


Ql  = 


-PL  flg- 

1  -  Qsa,  Dt 


if<7  <qsa,orDqsal  <0 
if4^M,orD<7M,  >0 


(3.18) 


with  L  the  latent  heat  of  vaporization,  q  is  the  water  vapor  mixing  ratio,  and  qsat  is 
the  saturation  mixing  ratio.  Equation  (3.18)  can  be  approximated  by 
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Ql  ~ - LpDc i  /  Dt  since  qsat  is  small  compared  to  unity,  and  Ql  =  0  when  Dq/Dt  =  0 
(Gill  1982).  The  heat  arising  from  frictional  dissipation  is: 


Q,  =  pv 


dx 


V  J  J 


(3.19) 


where  v  is  the  kinematic  viscosity,  and  the  indices  i  and  j  represent  summation  over  the 
three  spatial  dimensions.  Equation  (3.17)  then  becomes: 


DT  Dq  Dp 

pcD +  pL —  =  —  +  pv 

p  Dt  Dt  Dt 


\dxJJ 


+  kV2T  +  QR 


(3.20) 


where  the  mass  continuity  equation  has  been  used  to  put  the  first  two  terms  into  flux 
form.  This  equation  describes  the  internal  (sensible  and  latent  heat)  energy  of  an  air 
parcel.  To  derive  the  total  energy  equation,  expressions  for  mechanical  (kinetic  and 
potential)  energy  are  also  needed  for  the  parcel.  Starting  with  the  momentum  equation  in 
vector  form: 


Du  -  _ 

—  +  2£lxu 
Dt 


=  -Vp  -  pVO  +  Fu 


(3.21) 


where  u  is  the  three-dimensional  wind  vector,  2Q  x  u  is  the  Coriolis  acceleration,  g  is 
gravity,  and  Fu  is  molecular  friction/dissipation.  Dotting  Ft  into  (3.21)  yields  the  kinetic 
energy  equation,  where  Fu  has  been  separated  into  its  diffusive  and  dissipative  parts: 


2  2  \ 

+  v  +  w  I 

Dt 


-Fi-Vp  -  up  VO  +  ppV 


(  l  3 
1  -»2 

U~ 

\-pv 

V  2  J 

(3.22) 


with  p  denoting  the  molecular  viscosity.  The  potential  energy  is  given  by  the 
geopotential: 


DO  _  - 

p - =  pw-VO  =  wgp 

Dt 


(3.23) 
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Since  <J>  depends  only  on  z,  the  partial  derivative  with  respect  to  time  is  zero.  Adding 
(3.23)  to  (3.22)  yields  the  mechanical  energy  equation: 


1  D(if  +  v2+w2) 

-p— - - 

2  Dt 


D<b 


=  p/jV 


fi_0 

-  pv 

2 

+  | 

'  dp  Dp  > 

u  J 

v  dX  J  V 

v  3/  Dt  j 

(3.24) 


where  the  pressure  advection  has  been  replaced  using  the  definition  of  the  material 
derivative,  Dp/Dt  =  dp/dt +  u-Vp .  Adding  the  internal  energy  equation  (3.20)  and  the 
mechanical  energy  equation  (3.24)  yields  the  total  energy  equation  in  material  form: 


DT  T  Dq  1  D(u-+v2  +  w2) 

CD - +  Lp-  +  -p - 

'  Dt  Dt  2  Dt 


+  P 


Dm 

Dt 


_D 

Dt 


( 


c  T  +  Lq  H — "(m  +  v"  +  w~ )  +  O  —  ppV 


(im 

—u 

v2  j 


+  kS/'T  +  Qr  + 


(3.25) 


Note  that  the  dissipative  heating  terms  have  cancelled  since  they  represent  a  conversion 
between  mechanical  and  heat  energy.  Here,  it  is  assumed  that  the  diffusive  component  of 
friction,  thermal  conductivity,  and  radiation  on  the  right  side  of  (3.25)  are  small  and  can 
be  neglected  in  this  application.  The  time  derivative  of  p  on  the  right-hand  side  does  not 
cancel  because  the  thermodynamic  equation  was  expressed  in  terms  of  the  enthalpy 
instead  of  the  internal  energy  (e.g.,  cpT  =  cT  +  p/p).  In  this  application,  local  changes 

in  pressure  are  associated  with  acoustic  waves  in  the  atmosphere  and  this  term  is  small 
and  also  neglected.  With  the  right  side  set  to  zero,  Equation  (3.25)  is  known  as 
Bernoulli’s  equation  (Gill  1982).  Expanding  the  material  derivative  to  flux  form  in 
azimuthally-averaged  cylindrical  coordinates: 
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p 


c  T  +  Lq  +  —  (w2  +  v2  +  w2)  +  gz 
2  J 


dt 


pw 


pru 


cpT  +  Lq  +  —(u  +  v~  +  w  )  +  gz 


dz 


=  0 


(  1  A 

c  T  +  Lq  +  —  {u2  +  v2  +  w2 )  +  gz 
V _ 2 _ )_ 

rdr 


(3.26) 


Each  variable  is  then  separated  into  mean  and  perturbation  quantities  as  before, 
where  the  mean  is  both  an  azimuthal  and  temporal  average,  and  are  Reynolds-averaged 
using  the  Boussinesq  approximation.  For  simplification,  several  terms  may  be  grouped 
and  defined: 


.2  ,  ..2  ,  ,  2  > 


J  2  '  ° 

(3.27a) 

-  p{u'2  +  v'2  +  M/2 ) 

2 

(3.27b) 

F,t  =  pCpii'T' 

(3.27c) 

Fzr  =  pCpw'T' 

(3.27d) 

Frq  =  pLu'q 

(3.27e) 

Fzq  =  pLw'q' 

(3.27f) 

where  E  is  the  total  energy  of  the  mean  flow,  e  is  the  turbulent  kinetic  energy,  and  the  F 
terms  are  radial  and  vertical  fluxes  of  temperature  and  moisture  in  which  the  subscripts 
indicate  the  covariances.  The  Reynolds  averaging  yields: 
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d(pE)  |  d(pruE)  |  cXpwE)  |  d(rFrT)  |  3F;7.  |  d(rFrq)  >  OF 


Z9 


dt  rdr  dz  rdr  dz 

dp  de  d(rue)  d(we )  d(ru'e)  d(w'e ) 

3/  3? 


rdr 


dz 


rdr 


dz 


rdr 


dz 


d(rwT)  d{uTrz)  d(rv  Trd)  d(vx,e) 


rdr 


dz 


+  ■ 


rdr 


3z 


(3.28) 


Note  that  terms  involving  Tijdui  /3x;  cancel  due  to  conversion  of  mean  kinetic  energy  to 

turbulent  kinetic  energy  (e.g.,  see  Lindzen  1990,  90-92),  but  flux  gradients  involving 
d(Tijuj )  /  dxj  remain.  For  lack  of  an  established  term,  this  is  coined  the  “shear  flux,”  as  it 

represents  a  loss  of  kinetic  energy  through  the  interaction  of  turbulent  momentum  fluxes 
and  shearing  flow.  Grouping  terms  and  integrating  over  the  control  volume: 

J  2^  [puE  +  FrT  +  Frq  ]  lr2  dz  -  J  2 nrx  [ puE  +  FrT  +  Frq  ]  lr  dz 
+ J/  [PWE  +  Kt  +  E,  ]  l2  2  Ttrdr  -  J'2  [ pwE  +  FzT  +  Fzq  ]  I,  In  rdr  = 


n 

d  Zi  J} 


z\ 


dip_  -pE-e) 
dt 


Inrdrdz  +  J  2nr2 


-ue  -  u'e  +  wt,_  +  vt. 


re 


t  dz 
J 


-  J  _  2nrx  | -ue  -  u'e  +  wTrz  +  vTrd  J  lr]  dz  +  J  | -we  -  w'e  +  uTrz  +  vtz0  J  L2  2  nrdr 

-  J  | "-we  -  w'e  +  uxrz  +  vtz6  J  L  2  nrdr 


(3.29) 


Rearranging  (3.29)  and  dividing  by  2n,  the  radially-integrated  enthalpy  flux  at  z.i 


is: 


31 


J’2  lF*  ]  L,  rdr  =  j"  [Fzt  +  F,  ]  L,  rdr  = 

J/  r2  [puE  +  Frk  +  ue  +  u'e  -  wTrz  -  vT rg  ]L dz 
-  j"  "  rx  [  puE  +  Frk  +  ue  +  u  'e-wTrz-vTre]\ri  dz 


+ 


| _  pwE  +  Fzk  +we  +  w'e  -  uTrz  -  vxz6  I  rdr 


-J;  [pwE 

•J  J 


+  we  +  w  e  —  u Trz  —  vtz8 


L  rdr 
z\ 


d(pE  +  e) 
dt 


fdrdz 


(3.30) 


It  is  then  assumed  that  the  vertical  moist  enthalpy  flux  at  zi  can  be  represented  by  a  bulk 
formula: 


F#  =  FzT  +  Fzq  =  pcp  w'T'  +  pLw'q'  =  CKp  I  uh  I  (k*-k) 


(3.31) 


where  k'  is  the  saturation  moist  enthalpy  at  the  sea  surface  (k*  =  cpTSST  +  Lq*),  and  q*  is 

the  saturation  mixing  ratio  at  the  surface.  Neglecting  the  radial  variation  of  Ck  over  the 
control  volume  yields  the  bulk  enthalpy  exchange  coefficient: 


CK  =  (j/  r2  [ puE  +  Frk  +  ue  +  u'e  -  wTr:  -  vr,.e  ]  lr2  dz 
-  J  rt  [  puE  +  Frk  +  ue  +  u'e-  wtrz  -  vT,.e  J  lr  dz 


+  P\pxvE 

Jr,  L 


+  F.  +  we  +  w'e  —  ut,._  -  vt 


ze\K  rdr 


-  J  pwE  +  we  +  w'e  -  uTn  -  vr,g  I.  rdr 


+ff 

Jz .  Jr, 


d(pE  +  e) 
dt 


rdrdz  /  J  [p  I  uh  I  (k*  -  k)]  l„  rdr 


(3.32) 


Similar  assumptions  as  in  the  momentum  equation  are  made  to  group  unresolved  terms 
into  a  residual: 
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CK  =  (p2  r2  [puE] l2  dz  -  j'2  n  [puE] \ri  dz  +  J  2  [pwE] \Zi  rdr-j*  [pwE]l  rdr  +  FSHEAR  ^ 
| "  [p  I  uh  I  (k*  -  A:)]  I,  rdr  +  R 


(3.33) 


where  the  transport  of  energy  by  the  interaction  of  turbulence  and  the  mean  wind  through 
the  lower  boundary  FSHEAR  =  I* "  \_urn  +  vt_0  ]  I.  rdr  is  written  separately  since  it  may  be 

J  r{  1 

resolved  indirectly  via  the  results  from  the  drag  coefficient  calculation.  Alternatively,  one 
can  obtain  the  average  enthalpy  flux  by  dividing  (3.30)  by  the  integrated  radius.  The 
unresolved  residual  is 


R  =  [\;/2[_Frk  +  ue  +  u'e-w'crz-vTre\K  dz 
rx  [/y,  +  ue  +  u'e  -  wxrz  -  vxr6  ]  I,,  dz 
+1;[^  +  we  +  w'e  -  utr.  -  vT_e  ]K  rdr 


]\z  rdr+  f  2  P 

d(pE  +  e) 

\ 

rdrdz 

J  1  *  Z\  Jr, 

dt 

> 

J  [p  I  I  {k*  -  /;)]  \ZI  rdr 


(3.34) 


As  with  the  momentum  budget,  accuracy  in  the  retrieval  of  the  enthalpy  exchange 
coefficient  requires:  (i)  minimal  errors  in  the  specification  of  the  axisymmetric 
thermodynamic  and  kinematic  fields;  (ii)  strict  adherence  to  mass  continuity;  and  (iii)  a 
small  magnitude  of  the  unresolved  residual  terms.  The  validity  of  these  assumptions  is 
presented  in  the  following  section,  using  numerically  simulated  tropical  cyclones’ 
budgets. 

C.  PROOF  OF  CONCEPT:  VERIFICATION  USING  SIMULATED  DATA 

To  test  the  accuracy  of  the  methodology  described  in  the  previous  chapter, 
comprehensive  momentum  and  enthalpy  budget  analyses  using  idealized  numerical 
simulations  were  undertaken.  The  first  objective  is  to  retrieve  bulk  drag  and  enthalpy 
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exchange  coefficients  with  known  magnitudes  using  numerically  simulated  observations 
of  tropical  cyclones.  The  second  objective  (to  be  described  in  Chapter  II.D)  is  to  use  the 
simulated  data  sets  in  a  budget  error  analysis. 

Two  numerical  weather  models  were  used  to  test  the  budget  methodology.  The 
first  was  the  two-dimensional,  axisymmetric  hurricane  model  developed  by  Rotunno  and 
Emanuel  (1987,  hereafter  RE87);  and  the  second  was  the  three-dimensional  Advanced 
Research  Weather  and  Forecasting  Model  version  3.0. 1.1  (hereafter,  WRF).  Since  the 
main  goal  of  the  numerical  modeling  effort  was  to  produce  idealized,  but  still  physically 
realistic,  hurricane  simulations,  an  effort  was  made  to  closely  replicate  the  initial 
conditions  and  parameterizations  of  the  two  models.  The  initial  atmospheric  conditions 
used  for  both  models  was  the  thermodynamic  background  defined  by  the  Jordan  (1958) 
mean  tropical  Atlantic  sounding  that  was  modified  to  be  in  thermal  wind  balance  (Smith 
2006)  with  the  initial  tangential  winds  specified  from  an  analytic  mesoscale  vortex 
(RE87,  Equation  37).  The  sea-surface  temperature  was  fixed  in  the  RE87  model  at  26.5° 
C,  and  at  28.0°  C  in  the  WRF  simulation.  Both  simulations  used  warm  rain  microphysics, 
with  a  single  category  for  all  liquid  water  but  variable  terminal  fall  speed  depending  on 
the  mixing  ratio  in  RE87,  and  separate  categories  for  cloud  and  rain  water  in  the  WRF 
model  using  the  Kessler  (1969)  microphysics  scheme.  The  RE87  model  used  Newtonian 
cooling,  and  radiation  was  not  included  in  the  WRF  simulation.  In  RE87,  a  uniform  1  km 
horizontal  (radial)  and  250  m  vertical  grid  spacing  were  used  in  a  1500  km  by  25  km 
domain.  All  other  settings  are  similar  to  those  from  the  4X  run  of  Persing  and 
Montgomery  (2003),  except  the  surface  fluxes  as  described  below.  In  WRF,  a  four-way 
nested  domain  was  used  with  a  fine  mesh  containing  205  x  208  gridpoints  at  1  km 
horizontal  grid  spacing,  which  tripled  with  each  successive  mesh  to  a  5400  km  square 
coarse  domain  with  27  km  resolution.  Fifty  vertical  levels  were  unevenly  distributed  with 
the  highest  resolution  in  the  boundary  layer  using  the  default  WRF  sigma  levels.  The 
WRF  model  was  run  for  a  total  of  8  days,  and  reached  steady  state  around  150  hours  into 
the  simulation.  The  RE87  model  simulations  reached  steady  state  from  100-150  hours 
and  were  extended  to  12  days  due  to  the  lower  computational  requirements.  A  summary 
of  the  model  settings  is  given  in  Table  1. 
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In  addition  to  the  fundamental  geometrical  and  dynamical  differences  of  the  two 
models,  the  surface  flux  and  sub-grid  turbulence  parameterizations  use  different 
approaches,  which  provided  independent  storm  structures  and  magnitudes  of  the  surface 
exchange  for  testing  the  budget  methodology.  In  the  RE87  model,  sub-grid  turbulence  is 
parameterized  by  a  local,  first-order  closure  that  relates  turbulent  fluxes  to  resolved  grid- 
scale  gradients  by  an  eddy-viscosity  assumption.  Since  the  model  is  axisymmetric,  all 
three-dimensional  motions  are  therefore  parameterized  by  gradients  in  the  azimuthal 
mean  flow.  The  resolved  deformation  flow  and  a  prescribed  mixing  length  determine  the 
eddy  viscosity  following  Smagorinsky  (1963). 

Surface  fluxes  in  the  RE87  model  are  given  by  bulk  aerodynamic  formula,  such 
that  the  bulk  surface  exchange  coefficients  can  be  prescribed  exactly.  In  an  effort  to 
determine  whether  the  proposed  budget  methodology  can  distinguish  effectively  between 
different  magnitudes  of  the  drag  and  enthalpy  exchange  coefficients,  five  simulations 
with  different  coefficient  magnitudes  and  different  Ck/Cd  ratios  were  performed.  The 
“high”  magnitude  simulations  used  the  Deacon  formula,  as  originally  prescribed  in  RE87: 

CD  =  1 .1  x  1(T3  +  4  x  1CT5(w2  +  v2)1/2 1^/2 

with  the  corresponding  Ck  given  by  a  defined  ratio  of  either  1.0  (“highl”)  or  0.7 
(“high0.7”)  to  the  drag  coefficient.  The  “mid”  magnitude  simulations  used  the  same 
Deacon  formula  at  surface  wind  speeds  up  to  35  m  s'1,  but  the  exchange  coefficients  then 
were  capped  at  higher  wind  speeds.  Note  that  the  surface  fluxes  were  not  capped  above 
hurricane  force,  merely  the  wind-speed  dependence  of  the  exchange  coefficients.  Two 
“mid”  simulations  using  Ck/Cd  of  1.0  (“midi”)  and  0.7  (“mid0.7”)  were  performed.  A 
single  “low”  simulation  also  was  performed  in  which  the  bulk  exchange  coefficients  were 
derived  from  direct  aircraft  measurements  of  the  turbulent  fluxes  reported  in  Black  et  al. 
(2004),  Drennan  et  al.  (2004),  and  Zhang  et  al.  (2008)  (see  Chapter  II.B  for  a  review  of 
these  results).  For  the  “low”  simulation,  Ck  was  set  to  1.1x10"  for  all  wind  speeds,  and 
the  Deacon  formula  Cd  was  capped  at  1.9x10"  .  These  five  RE87  simulations  give  a  range 
of  coefficients  for  testing  the  sensitivity  of  the  budget  methodology. 
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The  WRF  model  has  several  turbulence  closures  available.  The  modified  Yonsei 
University  Scheme  (hereafter  YSU)  (Noh  et  al.  2003;  Hong  et  al.  2006)  was  selected  for 
this  study.  This  scheme  employs  a  non-local  boundary  layer  parameterization  in  which  a 
parabolic  eddy  viscosity  profile  is  prescribed  to  match  the  estimated  model  boundary 
layer  depth  (defined  in  the  model  as  the  top  of  a  well-mixed  layer  of  virtual  potential 
temperature).  Additional  vertical  fluxes  that  account  for  transport  by  large  eddies  and 
entrainment  at  the  top  of  the  PBL  are  also  included  in  this  formulation.  Above  the 
boundary  layer,  the  YSU  scheme  uses  a  local  eddy  viscosity  assumption  which  is 
somewhat  similar  to  the  RE87  model  except  it  is  formulated  in  three  dimensions.  Davis  et 
al.  (2008)  described  modifications  to  the  original  YSU  scheme  for  hurricane  simulations 
in  a  similar  manner  to  the  “mid”  RE87  runs  described  above,  such  that  the  drag  and 
enthalpy  exchange  coefficients  do  not  increase  beyond  35  m  s'1.  In  Davis  et  al.  (2008), 
the  bulk  exchange  coefficients  were  not  modified  directly,  but  rather  were  adjusted 

through  the  surface  roughness  by  setting  z0  =  lOexp^-lO/n!  ’)  with  lower  and  upper 

limits  of  1.25  x  10"  m  and  2.85  x  10"“  m,  respectively.  This  modified  surface  flux  option 
was  selected  for  this  study,  but  no  modifications  were  made  to  the  WRF  source  code.  The 
drag  and  enthalpy  exchange  coefficients  were  then  inferred  from  the  output  u*  and 
sensible  and  latent  heat  fluxes  at  each  gridpoint.  The  retrieved  WRF  fluxes  were  slightly 
different  from  those  described  in  Davis  et  al.  (2008),  because  examination  of  the  v3. 0.1.1 
code  shows  a  modified  roughness  formulation  that  saturates  at  28  m  s'1  instead  of  35 
m  s"1  (the  -10  in  the  exponent  is  equal  to  -9  in  this  version).  Although  the  WRF  model 
includes  an  optional  one-dimensional  ocean  mixed  layer  model  for  these  simulations,  it 
was  not  included  in  this  study. 

The  maximum  low-level  winds  for  the  simulations  are  shown  in  Figure  10.  The 
RE87  winds  are  reported  at  125  m  altitude,  and  the  WRF  winds  are  from  a  diagnostic  10 
m  value.  All  six  simulations  have  instantaneous  low-level  wind  speeds  exceeding  70 
m  s"1,  so  that  all  of  these  modeled  storms  are  at  category  five  intensity.  After  a  ~4-5  day 
intensification  period,  all  the  simulations  reach  a  relatively  steady  maximum  intensity 
ranging  from  75 — 120  m  s"1.  The  weakest  storm  was  the  “high0.7”  simulation,  and  the 
strongest  was  the  “midi.”  The  three-dimensional  WRF  simulation  and  the  “low”  RE87 
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simulation  were  the  two  slowest  developers,  with  both  storms  peaking  at  150  hours  near 
100  m  s”1  before  a  slight  weakening  period.  The  “0.7”  CK/CD  ratio  simulations  are  weaker 
than  the  corresponding  “1.0”  simulations  for  the  same  flux  formulation,  which  is  in 
accord  with  the  predictions  of  E-PI  theory.  However,  the  large  difference  between  the 
“midi”  and  “high0.7”  simulated  intensity  is  counter-intuitive  in  terms  of  E-PI  theory, 
which  depends  only  on  the  Ck/Cd  ratio  and  not  on  the  absolute  magnitudes  of  the 
exchange  coefficients.  The  RMW  is  ~5  km  smaller  in  the  “midi”  simulation  than  in  the 
“highl”  (not  shown),  which  is  consistent  with  a  higher  tangential  wind  for  a  similar 
absolute  angular  momentum  distribution,  but  a  more  detailed  comparison  of  the 
differences  in  the  model  simulations  is  beyond  the  scope  of  this  study.  While  some  of  the 
simulated  storms  are  significantly  stronger  than  observed  tropical  cyclones,  they  are 
sufficiently  realistic  for  testing  the  momentum  and  energy  budget  methodology  outlined 
in  Chapters  II.B  and  II.C,  with  the  caveat  that  the  exact  magnitudes  of  various  terms  may 
be  over-  or  under-estimated  compared  with  actual  storms. 

A  representative  sample  of  the  simulated  structures  from  the  RE87  simulations  re¬ 
analyzed  by  the  SAMURAI  software  developed  for  this  study  is  shown  in  Figure  11.  The 
SAMURAI  analysis  is  described  in  detail  in  Chapter  IV.  The  analysis  is  shown  here  to 
facilitate  a  better  comparison  of  the  structure  with  the  observed  storms  presented  in 
Chapter  V.  A  comparison  of  the  raw  model  results  with  the  SAMURAI  analysis  is 
discussed  in  Chapter  II.D.4  below.  The  6-hour  averaged  axisymmetric  tangential  wind 
(Fig.  11,  top)  at  168  hours  from  the  “low”  simulation  has  a  maximum  of  -110  m  s"1  at 
1  km  altitude  and  20  km  radius,  with  a  well-defined  secondary  circulation.  The  top  of  the 
inflow  layer  has  a  steep  downward  slope  inside  the  RMW  that  approximately  matches  the 
level  of  the  maximum  tangential  wind,  and  outside  the  RMW  the  top  is  at  a  relatively 
high  altitude  of  2  km.  Upward  motion  and  outflow  maxima  are  found  just  above  the 
maximum  tangential  wind.  The  secondary  circulation  is  aligned  closely  with  the  angular 
momentum  and  energy  contours  in  the  eyewall  region,  with  a  relative  energy  maximum 
(-351  kJ  kg”1)  near  the  surface  just  radially  inward  of  the  RMW.  The  model  vortex  has  a 
similar  kinematic  and  thermodynamic  structure  as  in  previous  RE87  simulations  (Persing 
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and  Montgomery  2003).  Although  some  variations  in  the  structure  are  simulated 
depending  on  the  flux  configuration,  all  runs  are  qualitatively  similar. 

To  obtain  the  WRF  axisymmetric  structure,  a  SAMURAI  analysis  was  performed 
using  circulation  centers  derived  using  the  method  described  in  Chapter  V.C.2.  Since  the 
simulation  was  on  an  /-plane  with  no  environmental  flow,  the  wind,  pressure,  and 
circulation  centers  were  nearly  coincident  and  practically  stationary  over  the  period. 
Hourly  analyses  were  created  using  the  15-minute  model  output  with  independent 
estimates  of  the  center  position  from  alternating  east-west  and  north-south  slices 
analogous  to  the  “figure-4”  flight  pattern  used  in  aircraft  reconnaissance.  No  noise  or 
center  position  perturbations  were  added,  which  thus  yields  the  “best”  analysis  possible 
for  these  simulated  fields.  The  intent  was  to  emulate  the  analysis  procedure  used  with  the 
CBLAST  dataset  and  provide  a  best-case  scenario  for  understanding  errors  associated 
with  neglecting  the  residual  budget  terms  discussed  in  Chapters  III. A  and  III.B,  as  well  as 
small  errors  introduced  in  the  analysis  process. 

An  example  of  the  model  output  used  in  the  analysis  is  shown  in  Figure  12,  with 
the  1  km  altitude  fine-mesh  horizontal  wind  vectors  and  speed  and  simulated  50  dBZ 
radar  reflectivity  contour  at  144  hours.  The  eyewall  is  apparent  as  a  nearly  symmetric 
15 -km  wide  annulus  of  high  reflectivity  co-located  with  a  region  of  strong  winds  over 
80  m  s'1  and  a  peak  wind  speed  exceeding  120  m  s'1  on  the  inner  edge  of  the  high-rain 
region.  Some  radial  outflow  is  also  evident  in  the  wind  vectors  at  this  inner  reflectivity 
edge.  Minimal  structural  changes  are  simulated  in  the  eyewall  region  over  the  24-hour 
period  (144-168  hours)  used  for  the  budget  analysis. 

The  axisymmetric  SAMURAI  analysis  at  144  hours  for  the  WRF  simulated 
structure  (Figure  13)  is  generally  similar  to  the  RE87  structure,  but  is  notably  different  in 
several  details.  For  example,  the  wind  maximum  in  the  WRF  simulation  is  found  at  lower 
altitudes  and  radially  inward  of  the  location  in  the  RE87  “low”  simulation.  Furthermore, 
the  sloping  inflow  reaches  only  1-km  height  in  the  outer  part  of  the  domain,  and  is  ~20 
m  s'1  stronger  than  the  RE87  inflow  near  the  surface,  and  is  topped  by  a  stronger  outflow 
region.  Mostly  congruent  secondary  circulation  and  angular  momentum  and  energy 

contours  are  also  found  in  this  simulation,  but  the  energy  values  are  higher.  It  is  noted 
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that  the  density  is  lower  in  the  WRF  simulation  than  in  the  RE87  simulation  (not  shown), 
such  that  mass-weighted  energy  is  not  as  different  as  Figures  13  and  14  might  suggest. 
However,  the  axisymmetric  kinematic  and  thermodynamic  structures  of  the  RE87  and 
WRF  simulations  are  clearly  distinct  from  a  budget  standpoint.  While  this  is  an 
interesting  research  result,  the  objective  of  this  section  is  to  test  the  budget  method  with 
reasonably  realistic,  balanced  terms,  and  not  to  establish  the  fidelity  of  the  simulations  to 
a  real  tropical  cyclone.  A  systematic  exploration  of  the  differences  between  the 
axisymmetric  and  3D  models  is  therefore  deferred  to  a  later  study. 

Bulk  exchange  coefficients  diagnosed  from  the  model  output  for  the  “highl,” 
“midi,”  and  “low”  RE87  experiments  and  WRF  are  shown  in  Figure  14.  The  RE87  code 
was  modified  to  output  all  of  the  subgrid  terms,  including  the  residual  terms  and  surface 
stress,  heat,  and  moisture  fluxes.  Since  these  surface  flux  terms  were  prescribed  by  the 
bulk  aerodynamic  formula,  the  diagnosed  exchange  coefficients  match  the  analytic  values 
very  closely,  as  expected.  For  the  WRF  simulations,  the  equivalent  drag  and  enthalpy 
exchange  coefficients  were  diagnosed  using  the  bulk  aerodynamic  formula  and  the  output 
friction  velocity  and  surface  heat  fluxes.  These  Ck  and  Cd  values  are  slightly  lower  than 
the  10  m  analytic  values  since  the  fluxes  were  calculated  at  the  lowest  sigma  level  (~40  m 
altitude).  For  the  WRF,  the  drag  coefficient  (Fig.  14a)  levels  off  at  28  m  s'1,  with  an 
increasing  Ck  (Fig.  14b)  near  1.9x10’  for  wind  speeds  around  60  m  s'  .  This  flux  profile 
corresponds  reasonably  well  to  the  “midi”  flux  configuration  used  in  the  RE87 
experiments.  All  of  these  magnitudes  are  within  the  range  of  expected  values  for  the 
exchange  coefficients  in  real  tropical  cyclone  boundary  layers. 

Since  all  of  the  subgrid  scale  terms  could  be  recorded  from  the  RE87  simulations, 
it  was  also  possible  to  test  the  proposed  retrieval  method  when  all  terms  in  the  budget 
were  included.  This  test  for  a  representative  sample  from  the  “highl,”  “midi,”  and  “low” 
simulations  is  shown  in  Figure  15.  The  exact  match  of  the  derived  and  prescribed 
coefficients  is  as  expected,  and  serves  as  a  validation  of  the  retrieval  method  when  all  the 
terms  are  known.  These  results  were  independent  of  the  control  volume  used  for  the 
retrieval.  Since  the  goal  was  to  simply  validate  the  retrieval  methodology,  an  equivalent 
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test  was  not  performed  with  the  WRF  model  due  to  the  added  complexity  of  calculating 
and  storing  the  Reynolds  averaged  fluxes  in  cylindrical  coordinates. 


Table  1.  Summary  of  numerical  simulation  configurations 


Model 

RE87 

WRF 

Horizontal  grid  spacing 

1  km  radial 

1  km  fine  to  27  km  coarse 

Vertical  grid  spacing 

250  m 

50  uneven  sigma  levels 

starting  at  40  m  spacing 

Microphysics 

Warm  rain 

Kessler 

Surface  Layer 

“highl,”  “high0.7,”  “midi,” 

YSU  similarity  with 

“mid0.7,”  “low”  bulk 

hurricane  modification 

aerodynamic  formulations 

(isftcflx=l) 

PBL 

1st  order,  local  closure 

YSU  scheme 

Lateral  Boundary 

Jordan  sounding 

Jordan  sounding 

conditions 

SST  (C) 

26.5 

28.0 

Radiation 

Newtonian  cooling 

None 
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Figure  10. 


Simulated  maximum  low-level  wind  for  the  control  runs  used  in  sensitivity 
testing.  Low-level  winds  are  maximum  instantaneous  10  m  winds  diagnosed 
every  15  minutes  from  the  WRF  simulation,  and  are  maximum  instantaneous 
axisymmetric  125  m  winds  every  six  hours  from  the  five  RE87  simulations 
(see  inset). 
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RE87  'low'  1 68  hour  tangential  wind  (color),  radial  wind  (contour),  and  secondary  circulation  (vector) 


Angular  momentum  (color),  total  energy  (contour),  and  secondary  circulation  (vector) 
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Figure  11.  RE87  “low”  simulation  axisymmetric  structure  at  144  hours.  Top  panel 
shows  tangential  wind  (color,  m  s’1),  radial  wind  (5  m  s’1  contours),  and 
secondary  circulation.  Bottom  panel  shows  angular  momentum  (color, 
106m2s'1),  total  energy  (350  kJ  kg  1  +  1  kJ  kg1  contours),  and  secondary 
circulation  (vector).  Solid  contours  indicate  positive  values,  and  dashed 
contours  indicate  negative  values. 
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144  hr  WRF  simulation  wind  and  reflectivity  at  1  km  altitude 


i  i  i  1  i  i  i  i  i  i  i  i  i  Li  i  ljti  l~*-[  ur  i  uri  isr  i  Lr_i  isr.i  wL.x  lA.  i  lA  i  i/i.  i  i/x'. 
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Figure  12.  Wind  vectors  (scale  vector  at  bottom  right)  and  speed  (color  scale,  m  s'1)  at  1 
km  altitude  and  simulated  radar  reflectivity  50  dBZ  contours  at  144  hours 
from  the  WRF  simulation. 
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WRF  144  hour  tangential  wind  (color),  radial  wind  (contour),  and  secondary  circulation  (vector) 
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Figure  13.  WRF  axisymmetric  structure  at  144  hours.  Contours  as  in  Figure  11. 
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Figure  14. 


Diagnosed  (a)  Cd  and  (b)  Ck  from  five  RE87  simulations  and  WRF 
simulation  (see  insets  for  color  symbols)  used  in  the  sensitivity  analysis. 
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Figure  15. 
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Retrieved  exchange  coefficients  from  selected  RE87  simulations  (see  inset) 
using  budget  methodology  with  all  terms  included. 


D.  ERROR  ANALYSIS 

1.  Impact  of  Unresolved  Budget  Terms 

The  baseline  validation  of  the  methodology  indicates  excellent  agreement 

between  the  derived  and  prescribed  exchange  coefficients  when  all  budget  terms  were 

included.  However,  several  of  the  terms  cannot  be  estimated  with  the  CBLAST  dataset. 

Although  the  actual  magnitude  of  the  unresolved  terms  in  nature  is  unknown,  first-order 

estimates  obtained  from  the  RE87  output  are  shown  in  Figure  16.  The  root  mean  square 

values  were  determined  from  1008  samples  of  the  6-hourly  averaged  simulated  values 

from  168-288  hours  in  the  RE87  integrations  with  different  flux  configurations.  Control 

volumes  were  varied  in  size  and  shape  from  18.5-44.5  km  radius,  and  from  the  lowest 

model  level  (125  m)  to  1.625-3.125  km  altitude  to  get  a  representative  sample  of  the 

various  quantities.  This  analysis  indicates  that  the  leading  term  in  the  budget  is  the  flux  of 

mean  quantities  through  the  outer  surface  of  the  control  volume,  and  the  second  leading 

term  is  the  mean  flux  through  the  top  of  the  volume — both  of  these  terms  are  known 

quantities.  Although  the  subtraction  of  several  large-magnitude  terms  to  obtain  a  small 
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residual  generally  is  problematic  in  finite  precision  mathematics,  in  this  case  the  desired 
surface  fluxes  are  the  third  leading  term,  so  there  is  some  confidence  in  the  ability  to 
retrieve  this  quantity  with  reasonable  accuracy.  Radial  eddy  fluxes  were  the  smallest 
terms  in  both  budgets,  which  suggests  that  they  can  be  neglected  with  minimal  error  in 
the  inner  core  region.  The  two  mean  fluxes  across  the  inner  and  bottom  interface  are 
known,  with  the  latter  being  a  smaller  term  that  mainly  satisfies  integrated  mass 
continuity.  Unresolved  terms  of  numerical  significance  are  the  volume-integrated 
tendency/storage  term,  the  “shear  flux”  term,  and  vertical  eddy  fluxes.  Each  of  these  is 
discussed  in  more  detail  below. 

Whereas  the  simulated  storms  were  relatively  steady  in  their  intensities  over 
multi-day  timescales,  the  wind  speed  timeseries  (Figure  10)  indicates  that  there  was  non- 
negligible  variability  on  hourly  or  more  frequent  timescales.  Since  the  tendency  term  is 
integrated  over  the  entire  control  volume,  it  has  the  potential  to  be  a  significant  term  in 
the  budget.  To  examine  the  variability  of  the  tendency  term  more  carefully,  budget 
retrievals  were  performed  using  both  instantaneous  and  time-averaged  model  output.  A 
finite-difference  calculation  for  the  momentum  and  energy  tendencies  on  the  original, 
staggered  grid  from  consecutive  model  timesteps  was  found  to  be  close  to  the 
recalculated  value  from  the  momentum  (3.11)  and  energy  (3.28)  equations  using 
unstaggered  model  output.  A  power  spectrum  of  the  intensity  time  series  has  a  spectral 
peak  at  the  25-minute  period  (not  shown),  which  is  consistent  with  the  period  of  a  local 
inertial  oscillation  (Shapiro  and  Franklin  1995).  However,  the  magnitude  of  the 
oscillation  may  be  larger  in  the  axisymmetric  simulation  than  in  observed  storms. 
Experiments  with  different  averaging  periods  appeared  to  damp  this  oscillation,  but 
introduced  additional  errors  in  the  storage  term.  Most  notably,  the  unstaggered,  time- 
averaged  output  did  not  explicitly  enforce  mass  continuity  with  the  numerical  precision 
of  the  staggered,  instantaneous  gridpoint  values.  This  discrepancy  proved  to  be  a 
significant  source  of  error,  with  fluctuations  in  the  retrieved  exchange  coefficient 
magnitudes  exceeding  500%.  This  unacceptable  level  of  uncertainty  could  be  reduced  by 
strictly  enforcing  mass  continuity  through  a  recalculation  of  vertical  velocity  from  the 
given  divergence  field,  but  this  adjustment  led  to  unbalanced  fields  that  required  non- 
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negligible  tendency  terms  to  close  the  budget.  In  this  sense,  a  larger  tendency  term  can  be 
alternately  viewed  as  an  artificial,  residual  storage  of  either  momentum  or  energy  inside 
the  volume  resulting  from  errors  in  the  resolved  fluxes  across  the  top  and  sides  of  the 
control  volume.  These  volume-integrated  tendency  terms  also  were  found  to  be  spatially 
coherent,  and  were  manifest  primarily  as  a  sensitivity  to  the  size  and  shape  of  the  control 
volume  if  the  terms  were  neglected.  A  more  detailed  examination  of  the  sensitivity  to  the 
control  volume  is  presented  in  Chapter  III.E,  but  this  order  of  magnitude  analysis 
indicates  that  the  estimated  error  in  the  exchange  coefficients  from  neglecting  this  term  is 
d20-100%  for  moderate  Cd  and  Ck  values  (e.g.,  between  1  and  2.5  xlO"3). 

The  “shear  flux”  term  represents  an  integrated  loss  of  kinetic  energy  through  the 
lower  boundary  resulting  from  the  interaction  of  the  mean  wind  shear  and  turbulent 

stress,  Fshear  =  I* "  [wTr,  +  vTzf)  ]  I.  rdr  .  The  shear  flux  term  is  -40-60%  of  the  magnitude 

of  the  surface  fluxes,  which  results  in  a  significant  low  bias  if  neglected.  However,  the 
surface  stress  can  be  estimated  by  the  bulk  aerodynamic  formula  using  the  drag 
coefficient  derived  from  the  momentum  budget,  and  the  mean  wind.  This  correction 
works  well  with  the  numerically  modeled  budgets  when  the  drag  coefficient  is  specified 
exactly,  but  in  the  real  data  the  errors  due  to  uncertainties  in  the  surface  stress  and  mean 
wind  causes  uncertainty.  Sensitivity  tests  indicate  that  a  ±50%  error  in  Cd  translates  to  a 
±20%  error  in  the  magnitude  of  Ck  from  this  term.  To  avoid  adding  too  much  noise  to  the 
Ck  estimates,  the  mean  Cd  derived  from  all  six  missions  was  used  in  this  study  to 
estimate  the  surface  stress  used  in  evaluating  the  “shear  flux”  for  Ck-  Including  errors  in 
the  mean  wind,  ±20%  appears  to  be  a  reasonable  estimate  of  the  Ck  error  introduced  by 
estimating  this  term. 

Unresolved  vertical  turbulent  fluxes  at  the  top  of  the  control  volume  are  typically 
the  same  sign  as  the  surface  fluxes  in  the  lower  troposphere,  which  represents  an 
unresolved  flux  of  momentum  into  the  volume,  and  a  flux  of  energy  out  of  the  volume. 
This  leads  to  an  underestimation  of  both  the  drag  and  enthalpy  exchange  coefficients  if 
neglected  (see  Equations  3.14  and  3.32).  In  the  RE87  simulations,  these  eddy  terms 
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showed  a  strong  dependence  on  the  altitude  of  the  volume  top,  and  were  10-30%  of  the 
surface  flux  magnitudes  for  the  1.625-3.125  km  tops  used  in  Figure  15. 

For  the  CBLAST  analysis,  the  choice  of  volume  top  is  not  straightforward. 
Traditional  formalism  for  the  planetary  boundary  layer  assumes  that  the  magnitude  of  the 
turbulent  fluxes  decreases  to  zero  at  the  top  of  the  PBL,  but  an  examination  of  the  vertical 
structure  of  turbulence  using  the  CBLAST  dataset  (Zhang  et  al.  2009)  suggests  a 
separation  of  the  TC  boundary  layer  height  defined  mechanically  and  thermodynamically 
(Smith  and  Montgomery  2009).  In  between  the  outer  rainbands  of  Hurricane  Isabel,  the 
potential  temperature  mixed  layer  depth  was  400  m,  while  the  momentum  and  humidity 
fluxes  decreased  to  zero  at  700  m,  and  the  inflow  layer  extended  to  ~1  km.  A  similar 
depth  of  the  inflow  layer  at  the  eyewall  was  reported  in  Montgomery  et  al.  (2006),  but  the 
virtual  potential  temperature  mixed  layer  was  only  150  m  deep.  Smith  et  al.  2008  have 
suggested  that  the  PBL  has  different  characteristics  in  the  main  updraft  region  versus  the 
inflow  region  outside  the  eyewall.  In  the  eyewall,  turbulent  eddies  may  be  transported 
upward  where  flow  erupts  abruptly  out  of  the  PBL.  Under  such  circumstances,  one  might 
expect  these  turbulent  contributions  to  increase  the  magnitude  of  the  unresolved  turbulent 
fluxes  at  the  top  of  the  control  volume.  Examining  a  range  of  relevant  altitudes  is  a 
practical  approach  to  addressing  this  uncertainty.  Sensitivity  tests  described  in  Chapter 
III.E  suggest  that  volume  tops  below  and  near  the  top  of  the  inflow  layer  are  the  most 
appropriate  for  the  budget  calculation. 
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Figure  16.  Root-mean  square  magnitudes  on  a  logarithmic  scale  of  the  (a)  momentum 
budget  terms  and  (b)  energy  budget  terms  from  the  RE87  simulations  with 
different  flux  configurations  (see  inset).  The  colors  of  the  budget  terms  along 
the  bottom  indicate  which  terms  may  be  calculated  or  estimated  from  the 
CBLAST  dataset  -  green  are  known  terms,  red  are  unknown  terms,  yellow  is 
an  indirectly  estimated  term,  and  blue  is  the  desired  surface  flux  term. 


2.  Sea-Surface  Temperature 

The  sea-surface  temperature  plays  an  important  role  in  determining  the  enthalpy 
exchange  coefficient  by  determining  the  saturation  enthalpy  at  the  sea  surface 
( k*  =  c  Tsst  +  Lq*)  and  therefore  must  be  measured  accurately.  Direct  measurements  of 

the  SST  were  made  during  the  Hurricane  Fabian  missions  by  AXBTs  released  by  the 
NOAA  aircraft.  A  comparison  of  the  AXBT  temperatures  and  TRMM  microwave  imager 
satellite  estimates  from  the  3-day  period  prior  to  the  first  mission  is  illustrated  in  Figure 
17  and  summarized  in  Table  2.  The  TMI-derived  sea-surface  temperatures  are  near  28.5° 
along  the  track  throughout  the  analysis  domain,  but  the  AXBT  temperatures  are 
consistently  1-2  degrees  lower  at  the  same  locations.  Since  decreases  of  the  SST  are 
known  to  occur  during  and  after  the  storm  passage,  the  splash  locations  of  the 
measurements  were  plotted  in  a  storm-relative  coordinate  system  (Figure  18).  The  AXBT 
measurements  were  near  the  eyewall  (red  annulus)  with  the  exception  of  AXBTs  #4  and 
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#8.  On  2  September,  the  two  AXBTs  on  the  right  side  of  the  eyewall  are  ~1  C  lower  than 
the  one  released  on  the  left  side.  Given  the  relatively  good  azimuthal  sampling  around  the 
eyewall,  an  average  value  of  27°  C  was  chosen  for  this  mission.  On  3  September,  all  of 
the  AXBTs  near  the  eyewall  agree  to  within  0.4°  C,  and  are  consistently  -1.2°  C  lower 
than  the  TMI  temperatures.  A  slightly  higher  value  of  27.5°  C  was  utilized  on  this  day.  It 
is  interesting  that  the  general  location  of  Fabian  on  3  September  coincides  with  that  of 
Isabel  on  13  September.  Since  the  SST  is  estimated  at  27.5°  C  on  both  of  these  days,  it 
appears  that  some  of  the  ocean  cooling  that  occurred  during  Fabian’s  passage  recovered 
in  the  10  days  between  the  two  storms.  Only  a  single  AXBT  was  available  on  4 
September,  which  was  adjusted  upward  slightly  to  28°  C,  given  the  warmer  SST  field 
derived  from  TMI  in  this  region. 

A  comprehensive  analysis  of  the  SST  for  the  potential  intensity  estimates  for 
Hurricane  Isabel  in  BM08  used  satellite,  buoy,  and  airborne  radiometer  measurements 
(Figure  19).  Comparisons  of  the  3-day  TMI  product  with  the  AVHHR  estimates  are 
good,  and  lend  confidence  to  the  prior  estimates.  Given  the  high  bias  of  the  TMI 
estimates  in  the  Fabian  case,  the  29°  C  on  14  September  used  in  BM08  was  adjusted 
downward  slightly  here  to  28.5°  C.  The  28.5°  and  27.5°  C  sea  surface  temperatures 
obtained  previously  for  12  and  13  September,  respectively,  were  not  modified.  It  is 
possible  that  these  estimates  are  a  little  high  for  these  calculations,  given  the  lower 
temperatures  found  underneath  Fabian’s  eyewall,  but  without  additional  evidence  they 
are  the  best  estimates  available  of  the  SSTs. 

SST  errors  do  not  affect  the  retrieved  enthalpy  flux,  but  do  affect  the 
corresponding  exchange  coefficient  via  the  surface  saturation  enthalpy.  The  percentage 
errors  in  Ck  for  0.5°  C  SST  error  increments  from  both  the  WRF  and  RE87  output  are 
shown  in  Figure  20.  The  two  models  agree  quite  well,  which  is  largely  expected  since  the 
error  is  not  dependent  on  storm  structure  but  only  on  the  prescribed  saturation  enthalpy 
value.  It  is  reasonable  to  assume  the  error  is  not  much  larger  than  1°  C,  since  all  of  the 
SST  measurements  are  within  a  few  degrees  of  each  other,  and  are  directly  measured  in 

the  Fabian  case.  This  would  suggest  an  error  of  no  more  than  20%  in  the  derived  Ck- 
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Table  2.  Table  of  SST  Observations  for  Hurricane  Fabian 


AXBT  -  TMI  AXBT/TMI 


AXBT# 

D/  HH:MM:SS 

AXBT  SST 

TMI  SST  Difference 

Average 

1 

2/  17:25:24 

26.64 

28.95 

-2.31 

27.795 

2 

2/  20:02:29 

27.61 

28.65 

-1.04 

28.13 

3 

2/  20:17:40 

26.88 

28.5 

-1.62 

27.69 

4 

2/  20:49:46 

26.72 

28.8 

-2.08 

27.76 

5 

3/  17:12:06 

27.73 

28.95 

-1.22 

28.34 

6 

3/  17:38:25 

27.4 

28.65 

-1.25 

28.025 

7 

3/  19:15:23 

27.34 

28.65 

-1.31 

27.995 

8 

3/  20:34:07 

27.64 

28.5 

-0.86 

28.07 

9 

4/  19:44:36 

27.49 

29.1 

-1.61 

28.295 

TMI  Sea  Surface  Temperature  3-day  average  ending  09/02/2003 
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Figure  17.  SST  derived  from  TRMM  Microwave  Imager  satellite  (average  SST  over  31 
August  to  2  September  in  color),  and  AXBT  data  released  into  Hurricane 
Fabian.  Track  of  Hurricane  Fabian  (dashed  best  track,  from  2  to  5  September) 
is  shown  for  reference. 
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Figure  18. 
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Storm-relative  AXBT  splash  locations.  AXBT  numbers  correspond  to  Table 
2.  Red  annulus  corresponds  to  approximate  eyewall  location  at  30  km  radius 
from  the  storm  center. 
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Figure  19.  SST  derived  from  AVHHR  satellite  (average  SST  over  4  to  10  September  in 
color),  and  NOAA  WP-3D  downward-pointing  radiometer  (thin  line,  from 
~18Z  13  and  14  September).  Tracks  of  Hurricanes  Fabian  (dashed  best  track, 
from  2  to  5  September)  and  Isabel  (dashed  best  track,  with  thick  white,  solid 
line  indicating  analysis  periods  from  16-23Z  on  12  to  14  September)  are 
shown  for  reference.  (From  Bell  and  Montgomery  2008) 


52 


Sensitivity  to  SST  Errors 


Figure  20.  Retrieved  Ck  sensitivity  to  errors  in  sea-surface  temperature  using  RE87  and 
WRF  model  runs.  Units  are  in  percent  for  the  Ck  error. 

3.  Circulation  Centers 

It  is  well  known  that  an  inaccurate  center  can  lead  to  errors  in  the  Fourier 
decomposition  of  a  storm’s  kinematic  and  thermodynamic  structure  in  cylindrical 
coordinates  (Marks  et  al.  1992;  Lee  and  Marks  2000;  Bell  2006).  While  a  transformation 
of  the  data  into  cylindrical  coordinates  does  not  change  a  field  (except  for  the  singular 
point  at  the  origin),  the  incomplete  representation  of  the  field  by  a  finite  series  that 
truncates  higher-order  harmonics  does  introduce  error.  This  limitation  is  particularly 
important  when  only  the  wavenumber  zero  mode  is  retained,  as  aliasing  to  higher 
wavenumbers  is  manifest  as  errors  in  the  derived  axisymmetric  winds,  energy,  and  mass 
fields.  For  scalar  fields,  center  errors  are  manifest  as  radial  displacements  of  the 
measured  quantity,  such  that  sufficient  averaging  can  largely  eliminate  random  errors. 
For  vector  fields,  a  center  displacement  in  the  radial  direction  has  the  same  effect  as 
merely  moving  the  observation  inward  or  outward  from  the  displaced  center.  Sensitivity 
tests  indicated  that  this  type  of  error  leads  to  minimal  errors  in  the  budget  retrieval. 
However,  displacements  of  the  center  in  the  transverse  direction  lead  to  a  re-partitioning 
of  the  winds  into  the  tangential  and  radial  component. 
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It  can  be  shown  from  geometrical  considerations  in  a  polar  coordinate  system  that 
the  re-projected  wind  components  are  given  by  u  =  ur  /  r  +  ve  /  r  and  v  =  vr  /  r  +  ue  /  r , 
where  e  denotes  the  transverse  error  in  the  center,  u  and  v  are  the  actual  radial  and 
tangential  wind  components,  r  denotes  the  radius  of  the  observation,  and  tildes  indicate 

erroneous  quantities  with  f  =  \le2  +  r1  .  Since  r  /  r  ~  1  for  typical  center  errors  away 
from  the  eye,  the  percentage  errors  are  approximately  given  by  the  second  terms  in  each 
equation.  Given  that  v  »  u,  it  is  clear  that  the  largest  errors  would  be  found  in  the  radial 
component  on  the  inner  side  of  a  budget  control  volume.  For  example,  a  5  km 
displacement  of  a  60  m  s'1  tangential  wind  at  20  km  radius  would  yield  a  15  ms'1  error  in 
radial  wind.  Weighting  this  error  by  mass,  radius,  and  energy  or  angular  momentum  and 
integrating  over  the  column/ring  yields  a  large  error  in  the  resolved  flux  across  the 
boundary  of  the  control  volume.  In  cylindrical  coordinates,  integration  errors  are 
compounded  as  the  depth  of  the  volume  increases,  which  lends  some  support  for  the  use 
of  shallower  control  volumes  that  still  exceed  the  depth  of  the  boundary  layer.  The  fact 
that  the  error  is  radially  dependent  means  also  that  an  artificial  gradient  of  the  radial  wind 
is  introduced,  which  produces  an  erroneous  residual  flux  into  the  volume  and  inaccurate 
radial  divergence.  Since  changes  in  the  divergence  affect  vertical  motion,  the  error  also 
would  be  spread  to  the  flux  across  the  top  of  the  control  volume.  All  of  these  error 
implications  add  sensitivity  to  the  choice  of  the  control  volume. 

Sensitivity  tests  with  the  RE87  model  confirm  the  large  errors  described  above 
when  applied  to  an  individual  model  output,  whether  instantaneous  or  time-averaged. 
Standard  deviations  of  the  percentage  error  for  both  drag  and  enthalpy  coefficients  were 
-25%  per  km  center  error,  such  that  center  errors  of  4  km  had  a  ±200%  variation  at  the 
95%  confidence  level.  It  is  partly  for  this  reason  that  individual  radial  penetrations  and 
dropsonde  sequences  were  not  used  to  calculate  the  budgets,  since  a  center  error  can  be 
significant  on  a  single  pass.  Fortunately,  if  the  center  displacements  are  predominately 
random,  then  averaging  of  multiple  passes  with  independently  derived  centers  reduces  the 
errors  in  the  winds  and  thermodynamics,  which  suggests  that  the  analytic  results  are 
really  a  worst-case  scenario.  However,  this  also  indicates  that  reducing  center  errors  is 
critical  to  obtaining  meaningful  results  from  the  budget. 
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Since  a  variety  of  ways  exist  to  define  the  center  of  a  tropical  cyclone,  there  is  no 
unique  tropical  cyclone  center,  but  there  may  be  an  optimal  center  to  minimize  aliasing 
associated  with  a  truncated  Fourier  representation.  The  simplest  definition  of  a  TC  center 
is  where  the  earth-relative  wind  goes  to  zero.  For  a  moving  vortex,  the  wind  center  then 
depends  on  the  translation  speed  and  is  not  ideal  for  most  quantitative  applications.  The 
most  commonly  used  center  in  TC  studies  is  the  6-hourly  “best  track”  center  published  by 
the  National  Hurricane  Center.  These  centers  are  derived  from  a  variety  of  sources  and 
are  accurate  to  approximately  10  km  for  strong  hurricanes.  While  some  of  these  centers 
can  be  highly  accurate  when  aircraft  reconnaissance  is  available  near  the  6-hourly 
interval,  the  low  temporal  resolution  of  the  fixes  smoothes  out  any  high  spatial  resolution 
movement  such  as  trochoidal  oscillations  (Nolan  et  al.  2001).  Although  the  best  track  is 
sufficiently  accurate  for  comparison  with  forecast  tracks,  which  typically  are  concerned 
with  errors  >100  km,  it  is  not  optimal  for  deducing  axisymmetric  structures.  In  previous 
studies  using  the  current  dataset  (Montgomery  et  al.  2006;  Bell  2006;  BM08),  the  centers 
were  determined  by  the  method  described  by  Willoughby  and  Chelmow  (1982),  in  which 
flight-level  pressure  and  wind  observations  are  used  to  find  the  streamfunction  minimum 
for  a  axisymmetric  vortex  in  gradient  wind  balance.  In  the  absence  of  multiple  local 
minima  that  may  arise  from  transient  mesovortices,  this  dynamic  center  is  accurate  to  3 
km  and  is  the  preferred  method  for  high-resolution  centers  when  only  flight-level  data  are 
available.  A  third  center- finding  method  proposed  by  Marks  et  al.  (1992)  used  Doppler 
radar  data  to  determine  the  circulation  center  that  maximizes  the  axisymmetric  tangential 
wind  at  the  RMW.  Since  this  maximum  is  an  inflection  point  in  the  radial  gradient  of 
tangential  wind,  it  also  maximizes  the  vorticity  and  therefore  the  circulation  inside  the 
RMW.  An  alternate  center  definition  using  both  the  gradients  of  angular  momentum  and 
energy  also  was  considered,  but  was  difficult  to  implement  given  the  need  for  accurate 
gradients  of  thermodynamic  information.  In  the  case  of  a  stationary,  circular,  vertically 
aligned  vortex  in  gradient  wind  balance,  the  wind,  dynamic,  and  circulation  centers 
would  be  identical.  As  a  TC  translates  and  its  structure  departs  from  axisymmetry  and 
vertical  alignment,  these  centers  begin  to  diverge. 
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For  the  current  study,  the  availability  of  edited,  dealiased  Doppler  radar  data  made 
the  circulation  center  an  attractive  choice.  To  determine  the  circulation  center,  a  two- 
dimensional  “simplex”  search  (Neldar  and  Mead  1965)  on  a  gridded  dual-Doppler  wind 
field  (Mohr  1988)  with  1.5  km  horizontal  grid  spacing  finds  the  point  that  maximizes  the 
average  value  in  an  annulus  around  a  specified  radius.  A  3-km  annulus  was  used  in  this 
study,  and  an  azimuthal  average  spanning  at  least  60  degrees  was  enforced  to  prevent 
small-scale  wind  maxima  from  dominating  the  average.  The  simplex  search  compares  the 
average  tangential  velocity  at  the  three  vertices  of  a  triangle  whose  centroid  is  an  initial 
guess  of  the  TC  center.  The  algorithm  uses  reflection,  contraction,  and  expansion  of  the 
triangle  to  find  the  maximum  tangential  velocity.  To  ensure  a  global  maximum  was 
found,  the  simplex  search  was  initialized  at  16  locations  over  15  radii  bracketing  the 
estimated  center  and  RMW,  respectively.  A  mean  center  was  deduced  from  the  different 
simplex  solutions  at  1-km  altitude  for  each  radial  penetration  of  the  airborne  radar. 
Cubic-spline  interpolation  of  the  individual  centers  over  time  was  then  used  to  create  a  1- 
second  track.  The  track  was  reviewed  subjectively  for  outliers,  and  some  minor 
adjustments  were  made  after  identifying  centers  that  may  have  been  biased  by  local  wind 
maxima.  The  accuracy  of  individually  derived  circulation  centers  is  ~3  km  due  to  errors 
in  the  gridded  dual  Doppler  analysis  and  spatial  averaging  used  in  the  simplex  search. 
However,  the  use  of  a  broader  swath  of  radar-derived  wind  data  helps  to  eliminate  centers 
associated  with  mesoscale  vortices  that  distort  the  center  fixes  from  the  Willoughby  and 
Chelmow  (1982)  method.  Most  importantly,  since  these  centers  maximize  the  symmetric 
tangential  wind  by  definition ,  they  also  minimize  any  artificial  high-order  harmonics  in 
the  Fourier  decomposition. 

An  example  of  the  center-finding  method  for  a  radial  penetration  of  Hurricane 
Isabel  at  1650  UTC  12  September  2003  is  shown  in  Figure  21.  A  composite  radar 
reflectivity  from  the  lower  fuselage  5  cm,  C-band  radar  (color)  reveals  the  nearly 
symmetric  eyewall  with  values  exceeding  40  dBZ  that  is  superimposed  on  a  swath  of 
winds  at  1-km  altitude  from  the  tail  Doppler  radar  on  the  southwest  side  of  the  vortex. 
The  white  typhoon  symbol  indicates  the  location  of  the  radar-derived  circulation  center, 
and  the  white  circle  indicates  the  estimated  RMW.  Note  that  the  vectors  are  nearly 
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tangential  to  the  white  RMW  circle,  and  that  the  RMW  circle  also  follows  the  curvature 
of  the  outer  edge  of  the  eyewall  reflectivity  well  in  this  region.  Discrepancies  between  the 
circle  and  reflectivity  elsewhere  are  likely  due  to  temporal  and  spatial  averaging  of  the 
higher  beamwidth  (e.g.,  coarser)  lower  fuselage  data.  The  NHC  best-track  center 
interpolated  to  this  time  (black  closed),  and  dynamic  center  from  the  flight-level  data 
(black  open)  are  also  shown  for  comparison.  The  best  track  center  in  this  case  is  very 
close  to  the  circulation  center,  with  the  dynamic  center  displaced  to  the  southwest  in  the 
low  reflectivity  region.  It  is  possible  that  the  dynamic  center  was  influenced  by  one  of  the 
large  mesovortices  present  in  the  eye  at  that  time  (Kossin  and  Schubert  2004;  BM08). 
This  comparison  indicates  that  the  radar-derived  circulation  center  is  the  optimal  one  for 
the  axisymmetric  budget  calculations. 

The  NHC  best-track  and  the  high-resolution  tracks  derived  from  the  simplex 
searches  are  shown  for  all  six  missions  in  Figure  22.  As  expected,  the  tracks  are  nearly 
coincident,  but  the  higher  resolution  motion  is  apparent  in  the  colored  tracks.  The  largest 
discrepancy  appears  to  be  on  13  September,  where  the  differences  are  likely  due  to  the 
different  center  definitions,  but  are  still  well  within  the  large  eye  and  RMW  on  this  day.  It 
is  reasonable  to  assume  that  no  systematic  biases  are  present  in  the  tracks,  such  that 
random  errors  in  the  center  would  be  largely  minimized  by  the  variational  analysis 
technique.  Sensitivity  tests  with  minor  variations  to  the  observed  tracks  were  consistent 
with  this  hypothesis,  and  did  not  indicate  any  significant  changes  in  the  retrieved 
exchange  coefficients  for  displacements  on  the  order  of  1-3  kilometers.  The  error 
remaining  in  the  exchange  coefficients  after  averaging  the  observations  from  multiple 
radial  penetrations  is  estimated  at  ±20%. 
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Figure  21. 


Isabel  09/12  16:50-16:58  UTC 


1  I  1  I  1  I  I  I  I  ■  I  <jbz 

12  16  20  24  28  32  36  40  44 


Radar  analysis  of  Flurricane  Isabel  from  16:50-16:58  UTC  12  September, 
showing  reflectivity  composite  from  lower  fuselage  radar  in  color  with  dual 
Doppler  horizontal  winds  at  1-km  altitude  (vectors).  Hurricane  symbols 
indicate  centers  derived  from  radar  (white),  NHC  best  track  (closed  black), 
and  aircraft  winds  (open  black).  White  circle  indicates  radius  of  maximum 
wind  from  radar-derived  center. 


58 


Figure  22.  NHC  best  tracks  (black  and  gray)  and  radar-derived  tracks  (color,  see  inset) 
used  in  this  study  for  Hurricanes  Fabian  and  Isabel. 

4.  Gridded  Analysis 

To  integrate  the  kinematic  and  thermodynamic  integrals  numerically  in  Equations 
(3.15)  and  (3.33),  the  data  must  be  available  at  specific  grid  points.  For  the  previous  TC 
momentum  and  energy  budget  studies,  this  was  done  by  a  manual  map  analysis  that  was 
then  interpolated  to  regular  intervals.  This  analysis  was  largely  superseded  by  the 
development  of  objective  analysis  techniques  that  created  a  gridded  data  field  by 
weighting  the  observations  by  their  distance  from  a  particular  gridpoint.  Modern  data 
assimilation  techniques  have  continued  to  improve  the  analysis  of  data,  by  including 
observational  error,  a  priori  background  estimates  of  the  atmospheric  state  and,  in  some 
cases,  time-dependent  probabilistic  information  on  the  background  errors. 

For  the  current  study,  a  variational  analysis  technique  was  developed  based  on  the 
work  of  Ooyama  (1987)  and  modem  data  assimilation  methods.  The  objective  of  this 
method  is  a  maximum  likelihood  estimate  of  the  atmospheric  state  for  a  given  set  of 
observations  and  error  estimates.  A  detailed  derivation  and  description  of  the  SAMURAI 
(Spline  Analysis  at  Mesoscale  Utilizing  Radar  and  Aircraft  Instrumentation)  is  given  in 
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Chapter  IV.  Although  the  statistical  correlation  between  the  resulting  SAMURAI 
analysis  and  the  CBLAST  observations  is  very  good  and  yields  high  confidence  in  the 
derived  TC  structures  used  in  the  momentum  and  energy  budgets,  errors  arising  from  the 
gridded  representation  of  the  observations  must  still  be  considered. 

Although  variational  analysis  is  superior  to  objective  or  manual  analysis  in  a 
probabilistic  sense,  shortcomings  still  exist  in  the  practical  implementation  due  to 
incomplete  knowledge  of  the  observational  and  background  errors.  In  a  variational 
analysis,  the  derived  structure  conforms  to  new  observations  where  available,  and  relaxes 
back  to  the  a  priori  background  state  where  no  new  information  is  provided.  Given  the 
relative  lack  of  gridded,  axisymmetric  analyses  of  mature  TC  structure,  the  error 
estimates  for  the  background  state  are  uncertain.  For  the  current  study,  a  Gaussian 
background  error  covariance  was  assumed  with  an  influence  length  scale  based  on  the 
observation  density.  The  background  errors  were  purposefully  set  high  in  order  for  the 
analysis  to  conform  primarily  to  the  observations  where  available,  but  this  has  the 
detrimental  side  effect  of  making  the  analysis  unconstrained  in  data-poor  regions. 
Similarly,  the  background  error  length  scale  acts  as  both  an  effective  distance  for  the 
observations  and  as  a  spatial  filter;  a  large  value  helps  spread  the  information  provided  by 
an  observation  but  also  tends  to  overly  smooth  the  analysis.  An  appropriate  balance  must 
therefore  be  struck  between  these  trade-offs.  This  balance  was  assessed  by  varying  the 
specified  length  scale  and  subjectively  examining  the  resulting  gridded  fields  of  the 
numerical  model  output  and  the  CBLAST  analyses. 

With  the  exception  of  the  transverse  streamfunction  x//,  the  majority  of  the  gridded 
fields  were  found  to  be  relatively  insensitive  to  the  specific  choice  of  the  error  length 
scale  other  than  the  level  of  detail.  Since  the  secondary  circulation  depends  on  the 
derivatives  of  this  streamfunction,  even  relatively  smooth  fields  of  if/  can  produce 
substantial  variations  in  u  and  w  in  regions  where  the  derivative  is  not  properly 
constrained  by  either  the  observations  or  the  background  field.  This  deficiency  was  found 
to  be  most  significant  in  the  lowest  100  meters  of  the  analysis  domain,  where  the  signal 
from  many  dropsondes  was  lost  and  no  Doppler  radar  information  was  available.  An 

example  of  the  sensitivity  of  the  radial  wind  at  the  surface  to  dropsonde  coverage  and 
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three  background  error  length  scales  is  shown  in  Figure  23.  In  the  main  eyewall  region, 
the  radial  wind  profiles  are  very  similar  and  provide  a  good  least- squares  fit  to  the 
observations.  In  the  data  gap  regions  near  the  eye  (10-20  km  radius)  and  outside  the 
eyewall  (>50  km  radius),  the  solutions  diverge  significantly.  The  most  radical  is  the  4x 
(e.g.,  400  m  vertical,  4  km  radial)  length  scale,  which  has  a  large  oscillation  in  the  radial 
wind  in  the  data-poor  region.  However,  this  was  eliminated  largely  by  using  the  6x  filter. 
The  differences  between  the  4,  6,  and  8x  radial  winds  are  pronounced  at  the  outer 
boundary,  with  nearly  a  40  m  s'1  difference  where  no  observations  are  available  to 
constrain  the  solution.  This  is  partially  because  the  absolute  errors  in  \ff  increase  with 
radius  as  the  magnitude  of  the  streamfunction  increases,  and  partially  due  to  the  lack  of 
an  ideal  vertical  boundary  condition  on  the  streamfunction  near  the  surface.  Since  y/  =  0 
at  the  surface  constrains  the  value  but  not  the  derivative,  setting  the  first  or  second 
derivative  to  zero  introduces  an  artificial  structure  in  the  radial  wind.  In  this  study,  a 
third-order  derivative  constraint  is  used  to  damp  oscillations  of  t//  near  the  surface,  but  the 
magnitude  of  this  constraint  is  constant  over  the  domain.  These  tests  indicated  that  the  6x 
background  error  length  scale  is  a  good  trade-off  between  too  much  smoothing  and  the 
data  density  constraints,  and  is  adequate  for  the  current  study.  More  research  on  the 
proper  boundary  condition  for  the  cubic  spline  analysis  of  the  streamfunction  is 
recommended  for  future  applications,  with  alternative  spline  coefficients,  a  variable 
derivative  constraint,  or  anisotropic  background  error  covariance  possibly  improving  the 
near-surface  solutions  in  data-poor  regions. 

An  additional  factor  that  improved  the  analysis  was  the  incorporation  of 
Doppler  radar  data.  In  BM08,  the  Barnes  objective  analysis  of  Hurricane  Isabel  used 
the  equivalent  of  a  ~10x  filter  to  damp  noise  and  spread  the  dropsonde  information 
across  data  gaps.  By  incorporating  radar  data,  the  dramatic  increase  in  both  azimuthal 
and  radial  data  density  for  the  kinematic  variables  adds  significant  value  to  the  analysis 
and  helps  reduce  the  amount  of  filtering  required.  Without  the  radar  data,  the  surface 
radial  wind  sensitivity  (Figure  23)  extends  through  a  much  deeper  layer.  Ultimately, 
higher  quality  and  quantity  of  observations  lead  to  better  analysis,  such  that  the  highest 
confidence  in  the  gridded  fields  was  in  the  eyewall  region  where  the  best  data  coverage 
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was  available.  The  budget  control  volumes  were  therefore  centered  on  the  regions  with 
both  the  best  observations  and  the  highest  wind  speeds. 

As  an  additional  test  of  the  errors  associated  with  the  gridded  analysis,  the  RE87 
model  output  was  treated  as  simulated  observations  and  analyzed  with  the  identical 
SAMURAI  configuration  used  for  the  real  observations.  A  comparison  of  the  “mid0.7” 
wind  field  at  168  hours  is  shown  in  Figure  24.  The  analyzed  and  the  raw  model  fields 
compare  well,  and  the  structure  is  very  similar  to  that  in  the  “low”  output  at  the  same 
simulation  time  (c.f.,  Figure  11).  Whereas  the  raw  model  and  analyzed  tangential  winds 
are  nearly  identical,  the  analyzed  radial  winds  are  noticeably  smoother  than  the  raw 
model  output,  especially  in  the  peak  magnitude  of  the  inflow  and  outflow  at  the  RMW. 
The  effect  of  the  boundaries  is  also  apparent,  with  a  broader  region  of  strong  outflow  at  4 
km  in  the  analyzed  field  and  a  bend  in  the  radial  wind  contours  near  the  surface.  Since  the 
lowest  model  was  125  m,  the  specification  of  winds  near  the  surface  is  arbitrary  in  this 
case  and  is  not  necessarily  representative  of  an  error.  In  general,  the  analysis  faithfully 
reproduces  the  kinematic  structure  with  a  0.9984  statistical  correlation,  a  -1%  RMS 
difference  of  the  variables  other  than  the  radial  and  vertical  wind.  There  is  a  -28%  RMS 
difference  for  the  radial  and  vertical  wind,  with  most  of  the  differences  occurring  near  the 
boundaries.  The  6x  filter  length  is  perhaps  too  heavy  for  a  “perfect”  data  distribution,  but 
is  a  fair  trade-off  with  the  CBFAST  data.  Given  the  uncertainties  described  above,  it  is 
estimated  that  errors  arising  from  the  gridded  analysis  procedure  contribute  less  than  a 
20%  error  in  the  derived  exchange  coefficients. 
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Figure  23.  Surface  radial  wind  sensitivity  to  background  error  length  scale  and  gaps  in 
the  dropsonde  data  on  13  September.  The  dropsonde  wind  observations 
below  50-m  altitude  on  13  September  are  marked  by  squares,  and  three 
analysis  curves  with  different  background  error  covariance  length  scales  are 
shown  in  color  (see  inset). 
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(a)  Original  RE87'mid0.7' output  at  168  hours 
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Figure  24.  RE87  raw  model  output  from  “mid0.7”  simulation  at  168  hours  (panel  a) 
compared  to  SAMURAI  analyzed  fields  (panel  b).  Tangential  winds  (m  s'1 
shaded,  scale  at  bottom),  and  radial  winds  (5  m  s'1  contours)  with  dashed 
contours  indicating  inflow  and  solid  contours  indicating  outflow.  The  thick 
black  contour  indicates  zero  radial  wind. 

E.  DERIVED  EXCHANGE  COEFFICIENT  SENSITIVITY  USING 

SIMULATED  DATA 

The  previous  analysis  of  the  (i)  unresolved  budget  term  magnitudes  and  (ii) 
SAMURAI  gridded  analysis  method  suggest  that  errors  in  the  derived  bulk  exchange 
coefficients  could  be  significant  even  if  the  SST  and  circulation  centers  were  known 
exactly.  This  motivates  the  following  test  of  the  methodology  using  simulated 
observations  from  the  numerical  model  output  analyzed  by  SAMURAI  with  no 
instrument  or  center  errors.  Using  this  test,  the  following  questions  are  addressed: 

1.  How  precisely  can  one  estimate  the  exchange  coefficients? 

2.  Is  there  an  optimal  geometrical  configuration  (width  and  height)  for  the 
budget  control  volume? 
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To  more  accurately  estimate  the  minimum  uncertainty  in  the  coefficients  given 
“perfect”  observations,  a  series  of  sensitivity  tests  was  conducted  to  guide  the 
specification  of  the  optimal  configuration  of  the  budget  control  volume  when  real  data 
will  be  used.  The  6-hourly  averaged  output  from  168-288  hours  from  all  five  of  the  RE87 
runs,  and  the  hourly-averaged  “figure-4”  analysis  from  144-168  hours  of  the  WRF 
simulation  were  analyzed  with  the  SAMURAI  program  using  the  same  configuration  as 
done  for  the  real  data.  A  wide  variety  of  control  volumes  were  then  tested  and  compared 
with  the  known  exchange  coefficients. 

One  of  the  first  steps  was  to  determine  the  sensitivity  to  magnitude  of  the 
exchange  coefficient  using  the  re-analyzed  RE87  output  as  shown  in  Figure  25.  The  large 
quantity  of  model  output  and  extensive  variety  of  the  control  volume  parameter  space 
(e.g.,  variable  inner  radius,  width,  and  depth)  yielded  -21,500  retrievals.  From  this  test, 
generally  more  negative  percentage  errors  are  found  for  the  “high”  drag  coefficients 
(>3.5xl03),  since  only  a  few  samples  are  above  zero  (Fig.  25a).  A  slight  linear 
correlation  exists  in  the  largest  enthalpy  coefficients  (Fig.  25b),  which  generally 
corresponds  to  increasing  average  wind  speed  in  the  control  volume.  For  the  lower 
magnitude  exchange  coefficients,  no  particular  dependence  is  found  on  the  magnitude  of 
the  coefficient,  with  numerous  samples  containing  both  high  and  low  percentage  errors. 
Although  the  percentage  error  in  Figure  25  seems  to  be  a  useful  metric  for  examining  the 
problem,  absolute  errors  for  the  higher  magnitude  coefficients  may  be  on  the  same  order 
as  the  lower  values  of  the  exchange  coefficients.  Since  a  systematic  bias  clearly  exists  for 
the  “high”  exchange  coefficients,  these  were  removed  from  further  consideration  to 
ensure  the  magnitudes  correspond  more  closely  to  the  actual  retrieved  magnitudes  from 
the  CBFAST  dataset.  A  total  of  -12,900  retrievals  remained  after  the  “high”  simulations 
were  removed.  This  subset  of  the  data  does  not  have  a  systematic  low  bias,  and  has 
percentage  errors  that  correspond  to  absolute  errors  with  similar  magnitudes. 

The  sensitivity  of  the  drag  coefficient  to  the  control  volume  geometry  for  the 
“low”  and  “mid”  simulation  samples  is  shown  in  Figure  26.  The  sensitivity  to  the  control 
volume  depth  (Fig.  26a)  has  a  parabolic  trend  with  an  apparent  low  bias  for  almost  all 
depths.  However,  a  positive  peak  in  the  distribution  is  found  around  1500  m,  which  is  just 
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below  the  inflow  top  for  these  simulations.  The  sensitivity  to  the  control  volume  width 
has  a  decreasing  spread  of  the  uncertainty  with  wider  volumes,  but  also  has  a  consistently 
low  bias  at  the  widest  annulus  of  22  km.  A  useful  metric  for  combining  these  two 
parameters  is  shown  in  Fig.  26c,  where  the  width/depth  “aspect  ratio”  defines  the  shape 
of  the  volume.  The  spread  of  the  uncertainty  decreases  from  low  aspect  ratios  (e.g.,  taller, 
thinner  volumes)  to  high  (e.g.,  shorter,  wider  volumes),  but  also  indicates  a  consistently 
low  bias  at  the  higher  ratios.  The  sensitivity  to  the  location  of  the  inner  edge  of  the 
control  volume  (normalized  by  the  RMW)  has  a  general  low  bias  but  with  a  peak  around 
0.8  in  the  distribution. 

The  RE87  sensitivity  of  the  enthalpy  exchange  coefficient  (Figure  27)  has  similar 
patterns  with  respect  to  the  depth,  width,  and  aspect  ratio  of  the  control  volume  (Figures 
27a-c),  and  to  the  location  of  the  inner  radius  (Figure  27d).  The  parabolic  shape  of  the 
depth  sensitivity  distribution  peaks  at  a  slightly  higher  height,  and  is  more  symmetric 
about  the  zero  error  baseline.  The  other  three  panels  indicate  a  trend  toward  a  low  bias  at 
the  upper  end  of  the  tested  ranges,  but  not  as  large  of  a  low  bias  as  the  drag  coefficient 
distributions  at  the  same  values. 

The  control  volume  sensitivity  tests  for  the  drag  coefficient  (Figure  28)  and  the 
enthalpy  exchange  coefficient  (Figure  29)  with  the  WRF  simulation  have  similar  trends 
with  respect  to  increasing  depth,  width,  aspect  ratio,  and  inner  radius  as  in  the  RE87 
results,  but  with  different  bias  characteristics.  A  parabolic  shape  is  evident  in  the 
sensitivity  to  the  depth  of  the  control  volumes,  with  peaks  at  a  lower  altitude  near  1000  m 
and  400  m  for  Cd  (Figure  28a)  and  Ck  (Figure  29a),  respectively.  The  percentage  errors 
become  more  negative  for  increasing  volume  width  (Figures  29a  and  29b).  The  aspect 
ratio  sensitivity  (Figures  28c  and  29c)  appears  to  have  the  most  similarity  with  the  RE87 
results,  and  indicates  a  decreasing  spread  of  uncertainty  at  higher  ratios.  However,  a 
distinct  difference  exists  in  the  overall  bias  of  the  WRF  results,  with  most  of  the  Cd 
sensitivity  tests  having  positive  errors  and  the  Ck  sensitivity  tests  having  mostly  negative 
errors.  Since  the  WRF  model  had  two  sigma  layers  below  100  m,  an  additional  sensitivity 
test  was  performed  by  varying  the  lowest  level  of  the  control  volume  (e.g.,  zi  in 
Equations  (3.15)  and  (3.33)).  This  test  indicated  a  general  tendency  for  slightly  lower 
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exchange  coefficients  with  zi  equal  to  100  m  versus  at  the  surface.  However,  no 
discernable  high  or  low  bias  existed  for  either  level  (not  shown). 

One  of  the  most  significant  findings  from  this  series  of  tests  is  that  the  errors  were 
not  uniform  from  model  to  model.  No  single  control  volume  was  found  that  minimized 
both  the  bias  and  spread  of  uncertainty  for  both  the  drag  and  enthalpy  exchange 
coefficients  in  both  models.  It  is  probable  that  the  magnitude  of  the  individual  unresolved 
budget  terms  is  different  in  the  3D  WRF  simulations  versus  the  2D  RE87  simulations,  but 
the  overall  errors  were  comparable.  To  determine  whether  this  was  primarily  due  to  the 
differences  in  the  2D  versus  3D  simulations,  or  some  other  facet  of  the  simulations,  an 
additional  sensitivity  test  varying  the  vertical  mixing  length  in  the  RE87  model  was 
performed.  This  parameter  partially  controls  the  size  of  the  turbulent  eddies  in  the  RE87 
simulations,  and  was  set  to  200  m  in  the  five  original  simulations.  Lowering  this  value  to 
100  m  for  the  “low”  flux  configuration  produced  a  simulation  with  a  similar  overall 
intensity,  but  with  a  slightly  shallower  inflow  layer.  The  drag  coefficient  sensitivity  to 
control  volume  depth  is  compared  for  these  two  simulations  and  the  WRF  model  in 
Figure  30.  Although  the  two  RE87  simulations  have  similar  error  profiles,  a  shift  to  a 
shallower  peak  in  the  parabolic  profile  is  found  for  the  100  m  vertical  mixing  length 
simulation.  This  suggests  that  at  least  some  of  the  error  differences  are  due  to  differences 
in  the  simulated  structure,  and  not  solely  the  2D  versus  3D  geometry  of  the  simulations. 

A  common  feature  of  all  the  simulations  is  a  trend  toward  a  decreased  spread  of 
the  percentage  error  for  larger  aspect  ratios.  This  can  be  interpreted  as  an  increased 
numerical  stability  as  a  larger  integrated  surface  flux  term  is  solved  with  wider  control 
volumes.  A  trade-off  then  exists  since  wider  volumes  must  account  for  the  variability  of 
the  wind  speed  with  radius,  but  given  the  larger  numerical  sensitivity  of  the  exchange 
coefficients  retrieved  from  thinner  volumes,  these  wider  volumes  appear  to  be  necessary. 
The  sensitivity  tests  also  indicate  a  parabolic  error  trend  for  increasing  depth  of  the 
control  volume,  with  a  peak  in  the  distribution  just  below  the  top  of  the  inflow  layer. 
Although  the  peaks  were  not  necessarily  correlated  with  a  zero  bias,  they  are  believed  to 
be  related  to  the  levels  at  which  the  vertical  eddy  fluxes  become  small  (i.e.,  the  top  of  the 
boundary  layer).  These  vertical  eddy  fluxes  may  not  go  to  zero  like  in  the  rainband  region 
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as  studied  by  Zhang  et  al.  (2009).  However,  it  is  reasonable  to  assume  that  the  vertical 
eddy  fluxes  are  smaller  at  or  near  the  top  of  the  boundary  layer  compared  to  the  surface 
fluxes.  Sensitivity  analysis  of  the  retrieved  fluxes  from  the  observations  show  similar 
trends  in  the  magnitude  of  the  fluxes  with  depth  of  the  control  volume.  The  parabolic 
profile,  as  opposed  to  an  increase  with  depth  followed  by  a  leveling  off,  suggests  an 
accumulation  of  errors  in  the  integrated  energy  or  momentum  flux  through  the  sides  of 
the  control  volume  with  increasing  height.  Therefore,  both  physical  and  numerical 
reasons  exist  to  suggest  that  control  volumes  with  tops  near  the  top  of  the 
boundary/inflow  layer  and  a  width  that  includes  a  sufficiently  large  integrated  surface 
flux  are  the  best  choice  for  minimizing  error. 

Since  the  volume  depth  is  constrained  by  the  inflow  layer,  the  aspect  ratio 
sensitivity  gives  reasonable  guidance  on  the  volume  width  necessary  to  reduce  the 
uncertainty  of  the  retrieved  Ck  and  Cd-  These  results  suggest  a  minimum  aspect  ratio  of 
~20  before  the  spread  of  the  percentage  error  decreases  to  ~±50%,  which  would 
correspond  to  a  width  of  20  km  for  a  1  -km  deep  volume.  The  data  distribution  ultimately 
constrains  the  maximum  width  of  the  control  volume  because  of  the  uncertainties  in  the 
structure  arising  from  fewer  observations  outside  of  the  eyewall  region.  Therefore,  an 
aspect  ratio  range  of  20-30  is  a  reasonable  range  of  the  control  volume  parameter  space 
given  the  observed  inflow  depths  and  data  distribution. 

A  summary  of  the  WRF  and  RE87  error  mean  and  standard  deviation  filtered  to 
these  aspect  ratios  is  shown  in  Figure  31.  The  absolute  error  (left  panel)  indicates  the 
smallest  bias  is  found  in  the  retrieved  WRF  Cd  and  RE87  Ck,  with  standard  devations  of 
0.76x10"  and  0.61x10’  ,  respectively.  A  distinct  low  bias  is  found  in  both  the  RE87  Cd 
and  WRF  Ck,  with  standard  deviations  of  the  error  of  0.49  and  1.04,  respectively.  In 
terms  of  percentage  error  (right  panel),  these  standard  deviations  correspond  to  21-55%, 
with  biases  ranging  from  +12  to  -59%.  The  largest  biases  or  standard  deviations  are  not 
confined  to  a  particular  model  or  exchange  coefficient,  but  appear  to  be  relatively 
independent.  Although  the  simulations  are  idealized  and  the  structures  are  only  semi- 
realistic,  it  is  believed  that  these  results  can  provide  useful  guidance  on  the  error 
characteristics  associated  with  a  particular  set  of  control  volumes,  but  do  not  specifically 
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describe  the  biases  of  the  real  dataset.  However,  the  significant  bias  and  uncertainty  of 
these  results  do  indicate  that  the  methodology  is  sensitive  to  the  control  volume 
specification,  even  within  this  range.  It  is  difficult  to  say  how  much  of  this  error  is  related 
to  the  models  versus  the  methodology.  It  is  reasonable  to  assume  that  some  of  this  error  is 
due  to  noise  in  the  models,  given  the  temporal  fluctuations  in  tangential  wind  apparent  in 
Figure  10.  However,  the  non-negligible  magnitude  of  the  unresolved  budget  terms  is 
likely  the  main  source  of  error.  The  errors  associated  with  the  control  volume  dimensions 
are  therefore  estimated  conservatively  at  ±50%,  but  it  is  acknowledged  that  they  may  be 
higher  than  this. 

A  summary  of  the  errors  identified  by  the  sensitivity  tests  is  given  in  Table  3.  The 
potential  errors  are  those  that  could  be  present  if  no  steps  were  taken  to  address  the 
problem,  with  the  estimated  errors  those  that  are  still  present  after  the  mitigation.  The 
largest  source  of  error  appears  to  be  the  neglect  of  the  unresolved  budget  terms, 
specifically  the  volume-integrated  tendency  and  vertical  eddy  fluxes  at  the  top  of  the 
control  volume.  This  error  is  estimated  at  ±50%  but,  as  mentioned  above,  this  may  be  an 
under-estimate.  The  remaining  errors  all  are  estimated  to  be  on  the  order  of  d20%  after 
mitigation.  Fortunately,  these  errors  are  independent  and  random,  and  there  is  no  reason 
to  expect  that  they  would  be  cumulative.  These  percentages  are  reasonable  estimates 
based  on  the  analytical  and  numerical  tests  performed  with  the  model  data.  These  tests 
also  have  identified  the  areas  for  which  additional  effort  was  required  to  ensure  the  best 
possible  results  from  the  observational  dataset.  If  left  unconstrained,  the  potential  errors 
listed  easily  overwhelm  the  magnitudes  of  the  retrieved  exchange  coefficients.  The  errors 
are  non-negligible  even  after  mitigation,  but  still  are  believed  to  allow  meaningful 
interpretation  of  the  results. 
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Figure  25. 
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150 


0.001  0.0015  0.002  0.0025  0.003  0.0035  0.004  0.0045  0.005 
Cr  Magnitude 


RE87  sensitivity  to  the  magnitude  of  the  exchange  coefficients.  Clusters 
indicate  similar  magnitude  coefficients  in  each  of  the  five  simulations.  Each 
dot  represents  the  percentage  error  of  an  individual  retrieval  of  an  exchange 
coefficient  using  a  single  control  volume  at  a  single  (6-hourly  or  hourly)  time 
interval  from  one  of  the  five  RE87  simulations. 
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Figure  26.  The  percentage  error  in  the  retrieved  Cd  from  the  RE87  simulations  obtained 
by  varying  (a)  the  depth  of  the  control  volume,  (b)  the  width  of  the  control 
volume,  (c)  the  aspect  ratio  (defined  as  the  width/depth),  and  (d)  the  location 
of  the  inner  radius  normalized  by  the  radius  of  maximum  winds.  Dots  are  the 
same  as  in  Figure  25. 
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Figure  27.  RE87  Ck  sensitivity  to  control  volume  specification.  Symbols  and  panels  as 
in  Figure  26,  except  for  Ck- 
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Figure  28. 


WRF  Cd  sensitivity  to  control  volume  specification.  Symbols  and  panels  as  in 
Figure  26,  except  for  WRF  Cd- 
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Figure  29. 


WRF  Ck  sensitivity  to  control  volume  specification.  Symbols  and  panels  as  in 
Figure  26,  except  for  WRF  Ck- 
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Figure  30.  Percentage  error  of  Cd  versus  control  volume  depth  for  RE87  “low” 
simulation  using  200  m  (green)  and  100  m  (red)  vertical  mixing  length  and 
WRF  (blue). 
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Figure  31.  (a)  Absolute  and  (b)  percentage  mean  bias  and  standard  deviation  for  control 

volume  aspect  ratios  of  20-30  for  combined  WRF  and  RE87  simulations  (see 
inset). 
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Table  3.  Summary  of  estimated  errors  in  budget  retrieval.  Error  percentages  are  valid  for  Cd 

values  from  1.9-2.4  xlO 3  and  Ck  values  from  1. 1-2.4  xlO"3. 


Error  Source 

Potential  Error 

Estimated 

Uncertainty 

Mitigation 

Mass  continuity 

residual 

>500% 

0 

Use  streamfunction 

for  analysis 

Shear  flux  term 

-60% 

±20% 

Estimate  term  using 

derived  Cd 

Other  unresolved 

budget  terms 

±200% 

±50% 

Constrain  control 

volume  size 

SST  Errors 

-20%  per  1°  C  error 

±20% 

Use  AXBT  data 

when  available 

Center  Errors 

±25%  per  km  error  for 

un-averaged  analysis 

±20% 

Use  high-resolution 

radar  circulation 

centers  and 

averaging 

Analysis  Errors 

Large  with  hand  or 

objective  analysis 

±20% 

SAMURAI  with  6x 

filtering 
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IV.  SAMURAI  ANALYSIS  TECHNIQUE 


Improvements  in  objective  analysis  techniques  for  gridding  meteorological 
observations  have  been  closely  coupled  to  the  increasing  volume  and  variety  of  data 
sources  used  in  NWP.  The  data  assimilation  problem,  in  which  the  best  initial  conditions 
for  a  weather  forecast  must  be  derived  from  sparse  observations,  is  essentially  identical  to 
the  analysis  problem,  in  which  the  best  estimate  of  the  atmospheric  structure  must  be 
derived  to  test  scientific  hypotheses.  In  both  cases,  the  researcher  or  forecaster  is 
attempting  to  determine  the  “true”  state  of  the  atmosphere  from  a  limited  set  of 
observations  that  contain  errors  that  result  from  sampling,  instrumentation, 
representation,  and  inversion.  An  analysis  that  is  used  to  produce  a  forecast  also  requires 
balance  enforcement  necessary  to  minimize  oscillations  of  the  primitive  equations,  versus 
a  stricter  adherence  to  the  observations  for  purely  diagnostic  studies. 

Early  objective  analysis  techniques,  such  as  those  proposed  by  Cressman  (1959) 
and  Barnes  (1973),  emphasized  the  primacy  of  the  observations  by  minimizing  the 
difference  between  a  gridpoint  value  and  nearby  observations  weighted  by  their  spatial 
distance.  These  well-tested  algorithms  still  are  commonly  employed  for  diagnostic 
studies  since  they  provide  a  simple,  functional  methodology  for  deriving  gridded 
atmospheric  structure  from  irregularly  spaced  data.  These  methods  have  been  largely 
superseded  within  the  NWP  community  by  variational  and  Kalman  filter  techniques, 
which  is  motivated  in  part  by  the  numerical  balance  requirements,  but  is  also  due  to 
advancements  in  the  application  of  probabilistic  theory  to  the  data  assimilation  problem. 
For  the  current  study,  these  advancements  are  adopted  within  a  mesoscale  analysis 
framework  in  order  to  deduce  the  most  probable  axisymmetric  TC  state  for  momentum 
and  energy  budget  calculations.  The  Spline  Analysis  at  Mesoscale  Utilizing  Radar  and 
Aircraft  Instrumentation  (SAMURAI)  software  program  has  been  developed  specifically 
for  this  study,  with  the  goal  of  obtaining  an  objective  analysis  that  can  incorporate 
aircraft  data  and  maximize  the  advantages  of  working  within  the  axisymmetric, 
cylindrical  geometry. 
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First,  a  variational  (VAR)  technique  is  derived  from  Bayesian  first  principles  to 
determine  the  maximum  likelihood  estimate  of  the  kinematic  and  thermodynamic 
structure.  While  3D  (or  2D  in  this  case)  VAR  approaches  are  not  new,  the  specific  design 
of  this  implementation  was  chosen  to  provide  an  optimal  framework  for  mesoscale 
analysis,  and  especially  for  tropical  cyclones.  Given  the  need  for  high-quality  vortex  and 
mesoscale  analysis,  particularly  for  initialization  of  numerical  models,  the  SAMURAI 
software  can  be  extended  naturally  into  asymmetric  cylindrical  and  Cartesian  geometries 
in  the  future. 


A. 


DERIVATION  OF  MAXIMUM  LIKELIHOOD  ESTIMATE 


The  derivation  starts  from  a  Bayesian  approach  following  Lorenc  (1986).  Bayes’ 
theorem  states  that  the  posterior  probability  of  a  state  A,  given  that  state  B  exists,  is 
proportional  to  the  prior  probability  of  state  A  multiplied  by  the  probability  of  state  B  if 
A  exists.  That  is: 


P(A  I  B)  oc  p(B  I  A)P(A) 


(4.1) 


where  state  A  is  defined  to  be  the  correct  estimate  of  the  atmospheric  state  given  some 
prior  observation  information.  In  this  case,  the  analysis  state  vector  x  is  equivalent  to  the 
true  state  of  the  atmosphere  xT.  The  prior  probability  that  this  estimate  is  true  is  then 

F(A)=  P(x  =  xt)  (42) 


Since  some  error  in  this  estimate  is  likely,  this  error  is  denoted  by  x£  =  xt  -  x ,  and 
the  spatial  covariance  of  these  “background  errors”  is  then  B  =  [x£xT£  ] .  Here,  the 

fundamental  assumption  is  that  the  errors  are  unbiased  and  the  probability  distribution  is 
Gaussian  about  the  background  estimate.  This  assumption  can  be  justified  by  the 
expectation  that  the  errors  are  random,  as  opposed  to  being  systematic.  The  Central  Limit 
Theorem  for  a  continuous,  multivariate  system  yields: 


P(x  =  xt)  = 


1 


ctn/Ztt 


exp 


f  o  -  M)2  A 

2cr2 


1 


IBU  yflK 


-exp 


(  i  T  x  \ 

--(.xt-xb)Tw  (xt-xb) 

V  2  7 


(4.3) 
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where  |B|  is  the  determinant  of  the  covariance  matrix  equivalent  to  the  standard  deviation 
a  in  the  univariate  distribution,  and  .17,  is  the  background  state  estimate  equivalent  to  the 
sample  mean  ju.  State  B  is  then  defined  to  be  the  state  where  the  set  of  observations  are 
obtained  by  sampling  the  (true)  atmospheric  state.  That  is,  P(B)  =  P(y  =  yt  I  x  =  xt) . 
Assuming  that  perfect  observations  could  be  obtained  by  some  forward  model  h(x)  that 
maps  from  the  true  state  to  observational  space,  then  the  error  is  denoted  by 
y£  =  yt-  h(xt  )  and  the  covariance  is  R  =  [y£yT£  ] .  These  errors  encompass  instrument  and 

representation  errors.  This  distribution  is  again  assumed  to  be  unbiased  and  Gaussian 
which  yields  a  similar  probability  distribution  function: 


P(y  =  y,  \x  =  xt)-- 


— 7 - exP\  -~(y,  ~  h(^t))TR  l(y,  -h(xt)) 

|R \iylln  V  2 


Given  these  two  probability  distributions,  Bayes’  rule  is  applied  to  derive  the 
probability  of  the  posterior  distribution  of  state  A  given  state  B;  e.g.,  the  probability  that 
the  atmospheric  state  is  correctly  estimated  given  that  measurements  of  the  true  state 
have  been  obtained.  Substituting  the  PDFs  into  Bayes’  rule  yields 

f  1  5 

P(x  =  x,  I  y  =  y, )  =  C  ■  exp  -  -  (x,  -  xb  f  B_1  (xt  -xh  )  +  (yl  -  h(x,  ))T  R_1  (yt  -  h(x, )) 

^  2  ^  (4.5) 


where  C  is  a  proportionality  constant  given  by 

P(B) 

1  1 

Rb  |b|2 


(4.6) 


Taking  the  natural  logarithm  of  both  sides  of  (4.5)  yields 


In  P(x  =  xt\y  =  yt)  =  \nC-  ~(xt  -  xb)  B"  (x,  -xb)  +  (y,  -  h(x,))  R”  (y,  -  h(x,)) 

V  l 


Thus,  the  probability  that  the  atmospheric  state  has  been  correctly  estimated  is 
maximized  when  the  variable  term  on  the  right  side  of  Equation  (4.7)  is  at  a  minimum. 
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Defining  the  term  in  parentheses  on  the  right  side  of  (4.7)  as  a  cost  function  J(x),  and 
dropping  the  subscript  t  for  clarity,  the  maximum  likelihood  estimate  (or  identically,  the 
minimum  variance)  is  achieved  when 

min{/(v)}=  min{/fc(x)  +  /0(x)}= 

min  (x  -  xb  f  B  1  (x  -  xb )  +  (y  -  h(x))T  R 1  (y  -  A(x))l 

L2  J  (4.8) 


This  minimization  is  achieved  when 

V/(x)  =  B-'(x  ~xh)  +  HrR_1(y -  h(x))  =  0  (4  9) 

where  Hy  is  the  linearized  Jacobian  matrix  of  h(x).  Note  that  if  h(x)  is  linear,  the 
Bayesian  solution  is  equivalent  to  the  least-squares  approach. 

The  form  of  the  cost  function  in  (4.8)  and  (4.9)  is  useful,  but  requires  the  storage 
and  inversion  of  the  background  error  covariance  matrix.  This  is  an  N  x  N  matrix,  where 
N  =  ixjxkxn  is  the  dimensional  space  of  the  analysis  given  by  the  product  of  the 
dimensions  of  the  spatial  vectors  and  number  of  variables.  This  inversion  can  be  avoided 
by  recasting  this  equation  in  an  incremental  form,  which  provides  an  improved  structure 
for  estimating  the  effects  of  this  matrix  without  requiring  its  full  storage  (Courtier  et  al. 
1994;  Huang  et  al.  2000).  Suppose  that  the  maximum  likelihood  estimate  can  be  obtained 
by  an  incremental  update  to  the  initial  background  estimate,  such  that 

x  =  xh+Sx  (4.10) 

Substituting  (4.10)  into  the  cost  function  (4.8): 

J(Sx)  =  -  ( Sx)T  B  1  (Sx)  +  (, h(Sx )  -  d)T  R"1  (h(x)  -  d) 

2  (4.11) 


where  d  =  y  —  h(xb)  is  termed  the  “innovation  vector.”  To  avoid  having  to  invert  B,  a 
symmetric  matrix  C  is  defined  that  has  the  same  eigenvectors  as  B  and  eigenvalues  of  the 

square  root  of  B,  such  that  B  =  Vfi  Vb  =  CrC  (Lorenc  1997).  A  control  variable  vector 
q  is  defined  as 
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q  =  C  l5x  =  C  1  (x  -  xb ) 

with  the  analysis  increment  in  physical  space  recovered  by 

8x  =  C  q 

Substituting  (4.13)  into  (4.11)  yields 

J(q)  =  ^qTq  +  ^  (HC^  -  d)T  R  1  (HC^  -  d) 

in  which  the  gradient  is  derived  by  differentiating  (4.14)  with  respect  to  q\ 

VJ(q)  =  (I  +  CrHrR“'HC)g  -  CrHrR“'J 


(4.12) 


(4.13) 


(4.14) 


(4.15) 


This  form  not  only  avoids  the  inversion  of  B.  but  has  several  other  benefits  as 
well.  The  identity  matrix  in  the  first  term  of  the  cost  function  prevents  the  smallest 
eigenvalue  of  the  Hessian  matrix  V2/  from  becoming  less  than  unity.  This  form  reduces 
the  condition  number  of  the  Hessian,  and  suggests  the  conjugate  gradient  algorithm 
should  converge  faster  than  the  form  in  (4.8)  (Gao  et  al.  2004).  More  importantly,  the 
control  variable  transform  allows  for  the  use  of  alternate  internal  state  vectors  in  different 
coordinate  systems.  Since  the  maximum  likelihood  estimate  is  highly  dependent  on  the 
specification  of  B,  a  form  is  desired  in  which  the  covariance  is  primarily  determined  by 
spatial  correlations  instead  of  cross-variable  correlations.  Thus,  variables  that  are  highly 
coupled  (i.e.,  u  and  w  in  the  tranverse  vortex  circulation)  can  be  solved  dependently 
through  the  use  of  a  single  variable  (i.e.,  streamfunction).  Likewise,  coordinate  or 
spectral  transforms  can  be  applied,  which  further  serves  to  isolate  the  background  errors 
by  effectively  diagonalizing  the  matrix.  Lastly,  since  the  background  error  covariance 
matrix  is  now  applied  in  a  forward  manner,  it  can  be  treated  as  a  linear  operator  that  can 
be  modeled  by  a  simplified  function,  such  as  a  Gaussian  filter.  This  avoids  the  storage 
requirements  of  the  N  x  N  matrix.  Specific  choices  related  to  C  used  in  the  SAMURAI 
technique  will  be  described  in  the  following  section. 
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B.  AXISYMMETRIC  SAMURAI  IMPLEMENTATION 


The  preceding  general  discussion  is  now  refined  for  use  with  axisymmetric 
hurricane  structure  and  available  CBLAST  observations.  Recall  that  the  probability 
distributions  associated  with  the  background  and  observation  errors  are  assumed  to  be 
Gaussian  and  unbiased,  and  the  inversion  of  B  has  been  avoided  by  use  of  the 
incremental  form  of  the  cost  function.  However,  some  description  of  the  C  and  the  h 
operators  and  the  control  variables  is  still  required.  The  background  error  correlations  are 
assumed  to  be  Gaussian  and  isotropic,  which  allows  for  an  efficient  recursive  filter 
operator  that  replicates  the  effects  of  this  correlation  (Hayden  and  Purser  1995).  For  the 
solution  to  be  the  most  effective,  control  variables  are  desired  in  which  the  analysis  fields 
are  smooth  and  well  approximated  by  this  spatial  error  correlation,  which  suggests  the 
use  of  the  following  control  variable  state  vector  in  cylindrical  space: 

q(r,z)  =  {prv,ll/,h',q'v,p'a}T  (4  16) 

where  p  is  the  total  density,  r  is  radius,  v  is  tangential  wind,  *//  is  the  transverse 
streamfunction,  h  is  the  saturated  moist  static  energy  given  by  CpT  +  Lqv  +  gz. ,  with  Cp 

the  dry  air  heat  capacity  at  constant  pressure,  T  the  temperature,  L  the  latent  heat  of 
condensation,  and  qv  the  water  vapor  mixing  ratio,  g  is  gravity,  pa  is  the  dry  air  density, 
and  primes  represent  departures  from  a  static,  background  reference  state.  All  of  the 
relevant  physical  variables  can  be  recovered  from  this  distribution,  and  these  quantities 
are  expected  to  vary  smoothly  on  the  chosen  grid.  Note  that  prv  is  chosen  instead  of  v  for 
two  reasons:  first,  to  match  the  physical  transform  of  the  streamfunction  into  mass- 
weighted  velocity  components,  simplifying  the  use  of  radar  data  by  a  uniform  weight  to 
the  radial  velocity;  and  second,  to  approximate  the  absolute  angular  momentum 
distribution  in  the  inner  core.  Thus,  the  first  three  variables  form  the  primary  set  for 
deducing  the  relevant  physical  quantities  for  use  in  the  budget  equations.  The  remaining 
two  quantities  are  used  to  isolate  the  density  from  the  mass-weighted  wind  components, 
and  temperature  from  the  moist  static  energy.  The  use  of  perturbation  quantities  allows 
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the  analysis  fields  to  avoid  strong  vertical  gradients  associated  with  the  background 
hydrostatic  pressure  and  water  vapor  distributions. 


To  convert  the  control  vector  q  from  gridded  variables  on  a  physical  grid  to 
analysis  increments,  the  background  error  covariance  is  formulated  as  an  operator 
sequence  of  a  cubic  spline  transform  that  is  followed  by  a  recursive  Gaussian  filter: 

B  =  crc  =  stbstastcscsasb  rztin 


where  Sb  is  a  transform  from  the  physical,  Gaussian  micro-mesh  (hereafter  referred  to  as 
the  “mish”)  to  a  nodal  representation,  Sa  is  a  transform  from  the  nodal  representation  to 
cubic  spline  coefficients,  and  Sc  is  the  recursive  Gaussian  filter  and  background  error 
inflation.  In  matrix  form, 

&  =  ScSASBq  =  DF(P  +  Q)'1  (A  m 


This  yields  a  set  of  spline  coefficients  in  physical  space  that  can  be  used  to  determine  the 
atmospheric  variables  and  their  spatial  derivatives  anywhere  in  the  domain.  Each  of  these 
operators  is  examined  in  more  detail  below. 

One  of  the  fundamental  aspects  that  distinguishes  the  SAMURAI  analysis 
technique  from  other  variational  solvers  is  the  use  of  a  Galerkin  approach,  which  is 
similar  to  the  Fourier  spectral  transform,  but  uses  the  cubic  B-spline  as  a  basis  (Ooyama 
2002).  This  has  the  advantages  that  the  interpolation  coefficients  for  a  function  or  its 
derivative  at  an  arbitrary  point  on  the  grid  can  be  easily  calculated  a  priori  for  each 
observation  or  balance  constraint,  instead  of  maintaining  a  memory  intensive  gridpoint 
representation  of  complex  interpolation  coefficients  and  finite  difference  operators.  Since 
the  basis  is  cubic,  it  is  computationally  efficient  and  continuously  differentiable  to  second 
order.  The  finite  element  representation  allows  for  a  flexible  incorporation  of  boundary 
conditions,  which  is  a  distinct  advantage  over  pure  spectral  techniques  for  non-periodic 
domains.  The  disadvantage  of  this  basis  is  that  it  is  not  orthogonal,  and  thus  requires  a 
matrix  inversion  to  obtain  the  spline  coefficients,  but  this  is  a  fair  trade-off  with  its  other 
desirable  characteristics. 
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The  spline  transform  (in  one  dimension)  is  given  by: 

a  =  (P  +  Q)-1b 

(4.19) 

where  a  is  the  state  vector  of  spline  coefficients  and 

u(r)=  Yj 

meM 

(4.20a) 

b  =  [bm]T ,  bm  =  jjm(rMr)dr 

(4.20b) 

P  =  ^mJr’  Pmm  =\D<t>»M)t(r)dr 

(4.20c) 

Q  =  f ,  qmm=\j:(rW:'{r)dr 

(4.20d) 

0n(r)  =  <b  ~ — —  >  for  m  g  M  and  r  e  D 
\  Ar  J 

(4.20e) 

Here,  u( r)  refers  to  a  continuous  representation  of  any  variable, 

Ct{r)  refers  to  the 

discrete  given  values  of  the  variable  (in  this  case,  the  control  variables  on  the  mish)  and 

O  refers  to  the  cubic  B-spline  given  by: 

hl-|«|  if  1*13*0 

o  J 

hz-13)1  if  2 > in > i 

u 

0  if  HI  >  2 

(4.21) 

The  one-dimensional  spline  and  its  derivatives  are  shown  in  Figure  32.  A  third 
derivative  constraint  is  incorporated  into  the  Sa  transform  (e.g.,  (P  +  Q)-1)  to  minimize 
Gibb’s  oscillations  and  spline  representation  errors  near  the  Nyquist  scale  (2Ax).  This 
also  acts  as  a  sharp,  6th  order  filter  when  used  repeatedly  as  in  the  conjugate  gradient 
minimization.  For  a  regular  r,z  grid,  the  inverse  of  (P  +  Q)-1  is  precomputed  via  a 
Cholesky  decomposition.  Since  the  spacing  of  the  internal  control  mish  is  arbitrary,  the 
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integrations  needed  for  b  and  (P  +  Q)  1  are  done  with  a  two-point  Gaussian  quadrature 
wherein  the  control  gridpoints  are  chosen  to  be  the  quadrature  points  of  the  evenly  spaced 
r,z  grid.  Note  that  the  Si  transform  (4.20a)  is  used  in  the  H  operator  converting  from 
spline  coefficients  to  observation  space.  Homogenous  or  inhomogeneous  boundary 
conditions  can  be  incorporated  easily  into  the  transform  by  a  direct  modification  of  the 
spline  coefficients  near  the  edge  of  the  domain.  For  this  study,  the  second  derivatives  of 
all  control  variables  were  set  to  zero  at  the  boundaries  except  for  \|/  and  prv.  The  variable 
prv  was  set  to  zero  at  the  axis  of  rotation  (r  =  0),  and  \|/  was  set  to  zero  at  the  axis  of 
rotation  and  the  surface  (z  =  0).  The  extension  to  two  dimensions  is  straightforward,  and 
more  details  on  the  cubic  spline  transform  and  the  effect  of  the  derivative  constraint  can 
be  found  in  Hausman  (2000)  and  Ooyama  (2002). 

The  operator  combination  Sc  =  DF  is  the  application  of  the  background  error 
covariance  matrix,  where  D  is  the  standard  deviation  of  the  background  error  and  F  is  the 
recursive  filter  given  by: 

n 

r>  =  Pp,  +  , 

i= i 

n 

^  =  A; +!>/,>, • 

;=1  (4.22) 

with  r  and  s  referred  to  as  the  advancing  and  backing  steps  of  the  filter.  It  can  be  shown 
that  this  filter  approximates  a  second-order  Gaussian  filter  with  minimal  computational 
cost.  This  formulation,  which  is  derived  in  Purser  et  al.  (2003),  improves  upon  previous 
multi-pass  filters  by  increasing  the  order  (n)  of  the  approximation  so  that  only  a  single 
forward  and  backward  pass  is  required.  Analytic  results  suggest  that  fourth  order  is 
sufficient  for  most  applications. 

The  resulting  state  vector  after  the  sequence  of  transforms  is  the  set  of  spline 
coefficients  for  the  state  vector  increments  on  the  physical  radius  grid.  Lastly,  the 
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variables  need  to  be  mapped  to  observational  space.  This  is  done  using  the  two- 
dimensional  form  of  the  cubic  spline: 

x(r,z)  =  11 

neN  meM  (4.23) 

Note  that  the  summation  occurs  only  over  a  few  nodes  near  the  physical  point  in  practice, 
due  to  the  finite  basis  representation.  The  basis  function  is  replaced  by  the  derivatives  of 
(f)m  for  the  streamfunction  to  convert  to  physical  variables: 


pru(r,z)  =  X  X  ~t,(r 

neN  meM 

prw(r,z)  -  X  X 

neN  meM 


(4.24) 


For  aircraft  and  dropsonde  observations,  the  direct  measurements  of 
thermodynamic  scalars  T,  p,  and  relative  humidity  are  transformed  to  h  ’,  qv,  and  p’a  in  the 
initialization  to  simplify  the  minimization.  Cartesian  wind  velocities  are  then  transformed 
into  their  cylindrical  projections  by: 

pu  =  p((w  -UJx  +  (v-  Vm  )y)  /  r 
pv  =  p((w  -  Um  )y-(v-  Vm  )x)  /  r 

pw  =  pw  (4.25) 

where  tilde  variables  are  Cartesian  wind  components,  v  and  y  are  Cartesian  distances  to 
the  storm  center,  and  r  is  the  radius  to  the  storm  center,  which  is  moving  at  (Um,Vm). 
Radar  data  requires  the  projection  of  the  analysis  variables  given  the  following 
transforms: 

pu  [(xsin(0)cos(y)  +  ycos(0)cos(7))  /  r]+ 
pv[vcos(0)cos(7)  -  ysin(0)cos(/))  /  r]+ 
pvv[sin(y)]=  pVD 

Vd=vd~  wt  sin(0)  -  Um  sin(0)cos(/)  -  Vm  cos(0)cos(7)  r4 


where  0  and  y  are  the  radar  beam  azimuth  and  elevation,  Vd  is  the  observed  Doppler 

velocity,  and  wt  is  the  estimated  terminal  fall  speed  of  the  precipitation  derived  from  the 
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reflectivity.  The  Joss  and  Waldvogel  (1971)  reflectivity /rain  relationship 
w™"  =  — 2.6 pZ 101  is  used  below  the  melting  level,  where  p  is  a  density  correction 
(Beard  1985)  and  Z  is  the  linear  reflectivity.  For  this  study,  the  top  of  the  analysis  domain 
was  restricted  to  2.5  km  such  that  the  rain  relationship  was  used  exclusively,  but  an 
ice/snow  relationship  (Atlas  et  al.  1973)  also  was  coded  for  use  in  other  applications. 

SAMURAI  was  coded  in  C++  with  the  cost  function  minimization  performed  by 
a  conjugate  gradient  algorithm  (Press  et  al.  2002).  Several  critical  loops  were 
parallelized  using  openMP  to  accelerate  convergence,  with  identical  results  to  serial 
execution.  For  the  current  study,  the  radial  nodal  resolution  was  set  to  1  km,  with  a 
vertical  resolution  of  100  m.  After  some  experimentation,  a  6Ax  length-scale  (6  km  radial 
and  600  m  vertical)  was  chosen  for  the  recursive  filter,  and  a  4Ax  scale  length  was  chosen 
for  the  third  derivative  constraint.  This  defined  the  approximate  spatial  influence  of  a 
single  observation  as  6Ax  and  damped  the  amplitude  of  spectral  features  less  than  4  times 
the  grid  spacing,  with  both  operators  effectively  acting  as  diffusive,  low-pass  filters. 

Single  observation  tests  using  synthetic  observations  were  performed  to  test  and 
debug  the  C++  code.  An  example  of  the  analysis  increment  using  a  single  Doppler 
radar  observation  is  shown  in  Figure  33.  The  instrumentation  error  of  this  test 
observation  was  set  very  low  to  induce  a  substantial  analysis  increment  for  illustrative 
purposes.  Since  the  Doppler  velocity  contains  information  about  all  three  components 
of  motion,  a  Gaussian  shaped  response  in  the  tangential  velocity  and  streamfunction 
results.  The  analysis  of  a  single  Doppler  velocity  observation  is  underdetermined  in  a 
traditional  objective  analysis  scheme,  but  in  the  variational  formulation  the  specified 
background  field  provides  the  additional  constraints  needed  to  produce  a  full  wind 
field.  The  use  of  the  streamfunction  as  the  control  variable  induces  a  toroidal 
circulation  around  the  observation  to  conserve  mass.  When  multiple  observations  are 
included,  the  result  is  the  maximum  likelihood  estimate  that  satisfies  the  Bayesian 
probability  derived  from  the  error  characteristics  of  the  observations  and  background 
analysis.  The  background  error  estimates  were  set  liberally  to  allow  strong  adherence  to 
the  observations  and  limited  dependence  on  the  background  field.  To  ensure  the  results 
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were  not  overly  sensitive  to  the  specified  background,  a  second  “outer  loop” 
minimization  was  performed  using  the  results  of  the  first  cost  function  minimization  as 
an  updated  background  field.  The  specified  background  errors  are  shown  in  Table  1. 
The  tangential  wind,  mass,  and  vapor  mixing  ratio  fields  from  the  Barnes  analysis  in 
BM08  were  used  for  the  Isabel  analysis,  with  simple,  analytic  fields  for  Fabian.  In  all 
cases,  the  transverse  streamfunction  background  field  was  set  to  zero. 

Though  the  error  characteristics  of  the  specific  aircraft  instrumentation  used  in 
this  study  are  well  known,  the  representation  errors  arising  from  the  use  of  an 
axisymmetric  coordinate  system  are  not.  Aircraft  observations  were  first  pre-processed 
into  a  moving,  storm-relative,  cylindrical  coordinate  system,  such  that  the  estimated 
observation  errors  were  a  combination  of  errors  in  both  the  instrumentation  and 
representativeness  in  2D-cylindrical  space.  A  full  description  of  the  instrumentation  and 
observation  errors  is  given  in  Chapter  V. 
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Figure  32.  (a)  The  cubic  B-spline  O  on  the  normalized  abscissa  c,  and  its  first  derivative 

O’;  (b)  the  second  and  third  derivatives.  From  Ooyama  (2002) 


Single  Observation  Test 
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Figure  33.  Single  radar  observation  analysis  increment  from  background  state.  Color 
indicates  increment  in  tangential  wind,  and  vectors  indicate  the  increment  in 
the  secondary  circulation. 
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V.  RESULTS 


A.  CBLAST  DATASET 

Six  intensive  observing  periods  (IOPs)  were  conduction  in  Hurricane  Fabian  from 
2-4  September,  and  in  Hurricane  Isabel  from  12-14  September,  2003  as  part  of  the 
CBLAST  and  NOAA/NESDIS  OCEAN  WINDS  experiments.  National  Hurricane  Center 
(NHC)  best-track  intensities  were  estimated  at  greater  than  120  kt  during  the  -16-23 
UTC  time  period  in  which  observations  were  collected  in  Fabian  (Figure  34)  and  Isabel 
(Figure  35),  which  makes  both  storms  major  hurricanes  on  the  Saffir-Simpson  scale 
during  the  six  IOPs.  Two  NOAA  WP-3Ds  (P3s),  the  NOAA  G-IV,  and  United  States  Air 
Force  (USAF)  Cl 30  aircraft  collected  in  situ  flight-level  and  dropwindsonde 
observations,  with  additional  Doppler  radar  and  radiometer  data  obtained  by  the  P3s  only. 
The  dataset  is  similar  to  that  used  in  Montgomery  et  al.  (2006),  Bell  (2006),  and  Bell  and 
Montgomery  (2008),  but  with  the  additional  use  of  the  University  of  Massachusetts 
stepped  frequency  microwave  radiometer  (SFMR)  and  Doppler  radar  in  the  axisymmetric 
composites,  and  the  addition  of  the  Fabian  observations. 

Flight-level  observations  at  one  Hz  resolution  from  the  NOAA  42  aircraft  used  in 
this  study  were  kindly  provided  by  the  NOAA  Hurricane  Research  Division.  A 
rudimentary  correction  for  instrument  wetting  errors  (Zipser  et  al.  1981;  Eastin  et  al. 
2002)  was  applied  to  supersaturated  dewpoint  temperature  measurements.  This  correction 
assumes  that  the  errors  for  the  temperature  and  humidity  sensors  are  equal  in  magnitude 
but  opposite  in  sign,  which  was  shown  by  Eastin  et  al.  (2002)  to  reduce  the  majority  of 
significant  wetting  errors  but  not  completely  remove  the  errors. 

National  Center  for  Atmospheric  Research  (NCAR)  Global  Positioning  System 
(GPS)  dropwindsondes  measure  pressure,  temperature,  relative  humidity  (PTH),  and 
horizontal  wind  speed  at  two  Hz  temporal  resolution  along  a  Fagrangian  trajectory  as  the 
sondes  fall  at  12-15  m  s’1  in  the  lower  troposphere.  This  yields  a  vertical  resolution  of 
approximately  5  m.  The  PTH  typical  errors  are  less  than  1.0  hPa,  0.2  C,  and  5% 
respectively,  and  the  wind  errors  are  less  than  2.0  m  s'1  (Hock  and  Franklin  1999). 
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Vertical  velocity  was  derived  by  removing  the  estimated  terminal  fallspeed  of  the  sondes 
as  a  function  of  pressure  (Wang  et  al.  2009).  All  dropsondes  were  quality-controlled  to 
remove  noise  and  other  instrument  errors  with  the  HRD  Editsonde  software.  Bell  (2006) 
provides  an  in-depth  discussion  of  the  dropsonde  data  from  the  CBLAST  experiment, 
including  issues  related  to  humidity  corrections  for  some  of  the  data. 

The  X-band  (3.2  cm)  wavelength  tail  Doppler  radar  employed  the  fore/aft 
scanning  technique  for  all  missions,  which  provided  radial  velocity  data  in  a  cone  ~20 
from  the  track  both  fore  and  aft  of  the  aircraft.  The  data  were  first  corrected  for 
navigational  errors  and  manually  edited  to  remove  ocean  returns,  radar  sidelobes,  and 
other  artifacts  (Oye  et  al.  1995;  Testud  et  al.  1995;  Bosart  et  al.  2002).  A  large  number  of 
radial  penetrations  were  edited  for  the  current  study,  which  included  all  of  the  multiple 
dropsonde  sequences.  The  times  of  the  edited  radar  data  used  in  the  analysis  are  shown  in 
Table  4.  Most  of  these  legs  were  used  for  circulation  center  fixes,  but  some  were  too 
short  to  obtain  meaningful  results  from  the  simplex  search  algorithm.  One  of  the  more 
challenging  aspects  of  utilizing  this  radar  dataset  was  the  velocity  ambiguity  resulting 
from  the  use  of  a  low  Nyquist  velocity  during  most  of  the  flights.  With  unambiguous 
radial  velocities  of  <20  m  s'1,  the  Doppler  information  frequently  was  “folded”  multiple 
times  in  conjunction  with  the  very  strong  winds  in  the  major  hurricane  eyewalls.  The 
Bargen  and  Brown  (1980)  algorithm  was  used  to  correctly  dealias  the  majority  of  the 
velocities  using  the  in  situ  aircraft  flight-level  wind  as  a  reference  point,  but  gaps  and 
noise  in  the  data  made  additional  manual  unfolding  of  many  rays  necessary.  An  iterative 
correction  was  made  for  poorly  dealiased  velocities  evident  in  the  dual-Doppler  wind 
fields  used  for  center  finding.  While  some  of  the  individual  radar  gates  still  may  contain 
dealiasing  errors  due  to  the  multiple  folds  in  high  gradient  areas,  they  are  not  likely  to 
affect  the  analysis  winds  significantly,  given  the  large  number  of  observations  and 
averaging. 

The  along-track  resolution  of  the  radar  scans  was  ~1.5  km,  with  a  75-m  range 
resolution  along  the  beam  and  1.9  degree  beamwidth.  In  an  effort  to:  (1)  provide  more 
uniform  spatial  resolution;  (2)  reduce  errors  associated  with  individual  gates;  and  (3)  thin 
the  voluminous  dataset  for  the  variational  analysis,  the  radial  velocities  were  averaged 
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along  each  beam  over  a  minimum  of  five  gates  that  increased  in  number  with  range  as  the 
beam  spread.  This  averaging  reduced  the  radial  resolution  to  approximately  the  diameter 
of  the  beam,  and  decreased  the  number  of  observations  by  nearly  an  order  of  magnitude. 
The  error  was  determined  from  two  contributions:  (1)  the  spectrum  width  of  the  radial 
velocity  measurements,  which  represents  the  variance  of  the  velocity  within  the  pulse 
volume  (Keeler  and  Ellis  2000);  and  (2)  the  estimated  error  in  the  terminal  fall  speed.  The 
sum  of  these  two  sources  with  a  minimum  value  of  2  m  s'1  was  used  for  the  error 
estimate. 

The  University  of  Massachusetts  deployed  a  SFMR  for  measuring  surface  wind 
speed  and  rain  rate  on  the  NOAA-42  aircraft.  This  C-band  (~5  cm)  wavelength, 
downward-pointing  radiometer  relates  brightness  temperatures  at  six  different 
frequencies  to  surface  emissivity  and  to  empirically  derived  wind  speeds.  Since  wind 
speed  is  a  nonlinear  observation  operator,  it  was  assumed  here  that  the  tangential  velocity 
comprised  the  majority  of  the  retrieved  wind  speed.  If  the  surface  inflow  angle  is  less 
than  30  degrees,  then  the  error  caused  by  this  assumption  is  no  more  than  -12%.  Given 
uncertainties  in  the  calibration  of  the  instrument  and  the  simplification  of  the  observation 
operator,  the  observation  error  was  set  to  a  relatively  high  value  of  10  m  s’1.  Nevertheless, 
these  SFMR  wind  speed  observations  provided  an  important  constraint  on  the  surface 
wind  speed  when  near-surface  dropsonde  winds  were  not  available. 

The  combined  instrument  and  representativeness  errors  for  the  observations  used 
in  this  study  are  shown  in  Table  5.  While  there  are  likely  some  larger  errors  associated 
with  individual  measurements,  these  values  are  believed  to  be  statistically  accurate  for 
the  bulk  of  the  observations.  The  analysis  composites  were  not  found  to  be  strongly 
sensitive  to  the  magnitude  of  the  prescribed  errors,  but  did  exhibit  a  general  trend  toward 
smoother  fields  as  the  instrument  error  values  were  increased.  If  the  errors  were  set  too 
small,  then  the  fields  exhibited  a  tendency  toward  over-fitting,  with  a  lower  root  mean 
square  difference  between  the  observations  and  analysis,  but  an  unrealistic  level  of  detail. 
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Figure  34. 


Date 


NOAA  Tropical  Prediction  Center/National  Hurricane  Center  (a)  best  track 
and  (b)  best  track  intensity  for  Hurricane  Fabian.  The  three  intensive 
observing  periods  on  2,  3,  and  4  September  are  highlighted  in  panel  (b). 
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Figure  35. 


Date 

NOAA  Tropical  Prediction  Center/National  Flurricane  Center  (a)  best  track 
and  (b)  best  track  intensity  for  Hurricane  Isabel.  Open  hurricane  symbol 
indicates  transition  to  tropical  storm  strength,  filled  symbol  indicates 
transition  to  hurricane  strength,  and  “L”  indicates  extratropical  transition.  The 
three  intensive  observing  periods  on  12,  13,  and  14  September  are  highlighted 
in  panel  (b). 
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Table  4.  Analysis  times  for  edited  Doppler  radar  data  for  Fabian  during  2-4  September  and 

Isabel  12  -  14  September  (all  times  UTC). 


Fabian  02 

Fabian  03 

Fabian  04 

Isabel  12 

Isabel  13 

Isabel  14 

1712-1726 

1704-1714 

1930-1949 

1650-1659 

1620-1642 

1618-1646 

1930-1941 

1839-1849 

2020-2035 

1721-1726 

1653-1709 

1716-1730 

1953-2002 

1849-1857 

2120-2139 

1828-1846 

1726-1742 

1816-1830 

2014-2027 

1929-1939 

2143-2159 

1901-1910 

1748-1807 

1849-1902 

2101-2114 

1940-1950 

2203-2219 

1943-1948 

1820-1841 

1906-1919 

2115-2123 

2010-2020 

2234-2250 

2006-2019 

1859-1913 

1936-1947 

2206-2219 

2100-2110 

2024-2036 

1916-1930 

1953-2009 

2223-2235 

2114-2124 

2038-2049 

1936-1950 

2013-2029 

2244-2252 

2204-2211 

2115-2130 

1953-2007 

2035-2050 

2026-2039 

2052-2111 

2043-2055 

2112-2125 

2102-2117 

2130-2147 

2151-2207 

2209-2225 

Table  5.  Combined  instrumentation  and  representativeness  errors  used  in  the  SAMURAI 

composites. 


Instrument 

(Type) 

pu  error 

pv  error 

pw  error 

q  ’  error 

p  ’  error 

h  ’  error 

Dropsonde 

(0) 

2  kg  m 2  s'1 

2  kg  nT2  s'1 

4  kg  m  2  s'1 

2  g  kg'1 

10  g  m 3 

5  kJ 

Flight  level 
(1) 

2  kg  m  2  s'1 

2  kg  m"2  s'1 

2  kg  m  2  s'1 

2  g  kg'1 

10  g/m"3 

5  kJ 

Doppler  Radar 
(2) 

Min(a+  2*sin(4>),  2.0)  kg  m'2  s’1 

SFMR 

(3) 

10  kg  m"2  s'1 
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B.  DERIVED  STORM  STRUCTURES 


SAMURAI  analyses  were  conducted  for  the  six  CBLAST  missions  using  the  data 
sources,  error  characteristics,  and  circulation  centers  described  above.  The  radar  and  non¬ 
radar  observation  data  distribution  for  the  2  September  mission  is  shown  in  Figure  36  in 
both  the  radius-height  and  polar  planes.  It  is  clear  that  the  radar  data  dominates  the  spatial 
distribution  due  to  the  wide  swath  of  the  tail  Doppler  radar.  Comprehensive  radar 
coverage  provides  good  confidence  in  the  kinematic  fields  throughout  the  domain,  except 
near  the  eye  where  scatterers  are  limited.  Excellent  dropsonde  coverage  is  also  obtained 
in  the  main  eyewall  region  near  the  RMW,  with  bracketing  observations  at  flight  level 
and  at  the  surface  from  the  SFMR.  Good  azimuthal  sampling  was  achieved  by  the 
multiple  sequences  across  the  eyewall.  Since  only  the  flight  level  and  dropsonde  data 
provide  thermodynamic  data,  the  highest  quality  energy  analysis  is  restricted  to  the  -20- 
30- km  annulus.  Note  that  the  closed  energy  contours  at  the  outer  radii  (>  -50  km  radius) 
are  likely  artifacts  from  the  lack  of  data  in  these  regions,  where  the  analysis  relaxes  back 
to  the  prescribed  background  state.  Similarly,  the  thermodynamic  structure  inside  of  20- 
km  radius  may  be  unreliable  due  to  a  data  void  region.  Control  volumes  were  not 
extended  to  these  regions  due  to  lack  of  confidence  in  the  thermodynamic  analysis 
outside  of  the  main  dropsonde  sequences.  However,  the  analysis  was  extended  to  the 
storm  center  and  just  beyond  the  main  dropsonde  sequence  locations  to  avoid  issues  with 
the  prescribed  boundary  conditions. 

A  large  number  of  observations  was  incorporated  into  the  analysis  composites.  A 
total  of  -711,050  observations  of  a  kinematic  or  thermodynamic  variable  were  included 
in  the  2  September  analysis,  with  the  largest  number  coming  from  the  radar  data.  A  bar 
chart  showing  the  relative  contributions  of  each  instrument  type  for  this  mission  is  shown 
in  Figure  37.  Nearly  600,000  Doppler  radar  velocities  were  included,  with  the  second 
largest  contribution  coming  from  dropsondes  at  -100,000.  Flight  level  and  SFMR 
observations  were  more  limited  due  to  their  single  altitude  sampling.  One  of  the 
advantages  of  the  variational  analysis  is  the  ability  to  effectively  combine  these 
observations  based  on  their  individual  error  characteristics.  The  statistical  comparison  of 
the  observations  versus  the  analysis  at  the  measurement  location  is  shown  in  Figure  38. 
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The  correlation  is  greater  than  0.99,  with  a  linear  fit  having  only  a  small  intercept  bias 
and  slope  near  1.  The  spread  of  the  difference  between  the  observations  and  analysis  is 
nearly  uniform  at  different  magnitudes.  The  units  of  the  observations  are  m  s"1  for 
kinematic  variables,  g  kg  1  for  water  vapor,  kJ  for  moist  static  energy  perturbation,  and  10 
g  m"  for  density.  Since  all  of  these  variables  have  similar  magnitudes  in  these  units,  the 
statistics  are  not  dominated  by  any  particular  measurement. 

The  prescribed  background  and  observation  errors  and  smoothing  parameters 
control  the  amount  of  spread  from  the  one-to-one  centerline.  The  distribution  of  the 
differences  is  shown  in  Figure  38b,  and  illustrate  an  approximately  normal  distribution 
centered  about  zero,  which  is  consistent  with  the  uniform  scatterplot  in  panel  a.  The  large 
number  of  observations  ensures  that  any  single  poor  measurement  will  not  bias  the 
analysis  significantly,  and  that  the  majority  of  the  differences  are  within  the  combined 
observational  error  and  variability  of  the  three-dimensional  observations  in  representing 
the  axisymmetric  mean.  These  statistics  suggest  a  good  fidelity  of  the  analysis  to  the 
observations  without  excessive  over-fitting. 

The  data  distributions  of  non-radar  observations  for  the  remainder  of  the  missions 
are  shown  in  Figures  39-43.  Radar  sampling  was  similar  on  all  days  to  the  first  mission 
and  is  not  shown.  The  dropsonde  distributions  are  similar  to  the  first  mission,  with  the 
exception  of  little  data  outside  of  40-km  radius  but  excellent  azimuthal  sampling  on  3 
September,  and  limited  azimuthal  sampling  on  4  September.  All  of  the  Isabel  missions 
are  characterized  by  comprehensive  radial  and  azimuthal  sampling.  Note  that 
distributions  are  slightly  different  from  those  shown  in  BM08  due  to  the  use  of  radar- 
derived  circulation  centers  instead  of  flight-level  winds  as  described  in  Chapter  III,  but 
the  patterns  are  qualitatively  the  same. 

The  derived  wind  fields  for  Hurricane  Fabian  are  shown  in  Figure  44.  The 
tangential  winds  depict  a  general  weakening  over  the  three  analysis  days,  with  a  decrease 
in  the  depth  of  the  strongest  winds  and  slow  weakening  of  the  maximum  tangential  wind 
by  the  third  day.  The  reduction  in  the  depth  and  intensity  of  the  tangential  wind  was 
concurrent  with  a  reduction  in  the  inflow  depth  and  magnitude  over  the  three  days.  The 

primary  updraft  also  appears  to  have  weakened  during  these  three  days.  The  RMW 
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remained  consistently  just  inside  30-km  radius,  with  a  sharp  gradient  of  tangential  wind 
towards  the  center.  Note  that  the  low-level  inflow  continues  past  the  RMW  on  all  three 
analysis  days,  with  the  flow  turning  upwards  near  the  high  tangential  wind  gradient 
region.  The  absolute  angular  momentum  and  total  energy  are  shown  in  Figure  45.  The 
secondary  circulation  generally  follows  the  contours  of  angular  momentum  above  the 
inflow  layer,  with  more  tilted  angular  momentum  surfaces  outside  the  RMW  on  the 
second  and  third  days.  The  energy  contours  also  depict  a  toroidal  shape  that  is  similar  to, 
but  not  congruent  with,  the  momentum  contours  and  secondary  circulation.  A  reduced 
radial  energy  gradient  is  apparent  on  each  consecutive  day,  which  is  consistent  with  the 
weakening  trend  in  the  kinematic  variables.  The  oppositely  directed  gradients  of 
momentum  and  energy  are  similar  to  those  found  in  the  numerical  simulations  in  Chapter 
III.E,  but  the  specific  magnitudes  and  structures  of  the  energy  are  notably  different  from 
the  RE87  and  WRF  simulations. 

The  analyzed  kinematic  fields  for  Hurricane  Isabel  for  all  three  missions  are 
shown  in  Figure  46.  Note  that  the  inner  radius  shown  is  15  km  for  clarity  due  to  the  larger 
domain,  although  the  analysis  domain  was  extended  to  the  TC  center.  The  qualitative 
features  of  the  SAMURAI  analysis  are  similar  to  the  Barnes  analysis  used  in  BM08,  but 
with  enforced  mass  continuity  of  the  secondary  circulation  and  a  notably  sharper 
resolution  of  the  boundary  layer  inflow.  The  characteristic  “over-shoot”  of  the  inflow 
past  the  RMW  is  apparent  during  all  three  missions  as  was  the  case  in  the  Fabian 
analysis,  with  a  sloping  inflow  top  that  coincides  with  the  height  of  the  maximum 
tangential  wind.  A  gradual  weakening  and  expansion  of  the  tangential  wind  as  described 
in  BM08  is  evident,  with  a  well-defined  outflow  and  updraft  core  just  above  the 
maximum  tangential  wind  on  all  three  days.  The  angular  momentum  and  energy  analyses 
are  shown  in  Figure  47.  The  consistent  structure  of  the  secondary  circulation,  momentum, 
and  energy  contours  is  similar  to  in  the  analyses  for  Fabian.  The  total  energy  structure  is 
also  qualitatively  similar  to  the  moist  entropy  structure  shown  by  BM08  (c.f.,  their  Figure 
5),  as  would  be  expected  by  their  similar  thermodynamic  variables. 

The  kinematic  and  thermodynamic  structures  apparent  in  the  six  analyses  are 
consistent  in  their  depiction  of  gradually  weakening  storms  just  past  their  peak  intensity. 

99 


Differences  exist  in  the  evolution  of  the  RMW,  with  Isabel  expanding  while  Fabian  remained 
nearly  constant,  and  in  their  respective  intensities.  These  six  missions  make  an  excellent 
dataset  for  the  flux  retrievals,  since  they  are  all  nearly  steady  state,  span  approximately  20  m 
s'1  in  their  peak  tangential  wind  speeds,  and  have  a  variety  of  different  inflow  depths  and 
magnitudes. 

The  analysis  statistics  for  all  six  missions  are  listed  in  Table  6.  The  number  of 
observations  varies  from  -500,000  to  more  than  1.7  million  on  14  September  (due  to  the 
large  domain  required  to  encompass  the  45  km  RMW),  with  comparable  differences  between 
the  analyses  and  observations  on  all  days.  Similar  minimum  and  maximum  differences,  low 
skewness,  and  kurtosis  values  near  3  suggest  nearly  normal  distributions.  One  exception  to 
this  statement  is  found  on  13  September,  where  the  high  kurtosis  and  larger  peak  differences 
are  most  likely  attributed  to  some  faulty  velocity  dealiasing,  but  the  low  mean,  RMS,  and 
standard  deviation  of  the  differences  suggest  these  are  isolated  errors.  Linear  correlations  are 
near  0.99  for  all  days,  with  low  intercept  and  slope  values  near  one  (not  shown)  for  all 
distributions.  These  statistics  suggest  that  the  SAMURAI  analyses  are  indeed  maximum 
likelihood  estimates  for  the  TC  structures  in  areas  of  sufficient  data  density.  Since  there  are 
no  high  confidence  a  priori  estimates  of  the  structure  available,  other  than  the  previous 
Barnes  analysis  with  the  same  dataset,  the  data  is  only  weakly  constrained  in  poor  data 
regions.  Fortunately,  the  rich  data  regions  encompass  the  areas  required  to  obtain  optimal 
aspect  ratios  of  20-30  for  the  control  volumes. 


100 


(a) 


(b) 


Fabian  09/02  Non-radar  observations 


Fabian  09/02  Radar  Observations 
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Figure  36.  Spatial  distribution  of  observations  on  2  September,  (a)  Non-radar  and  (b) 
radar  observations  in  the  radius-height  plane,  (c)  Non-radar  and  (d)  radar 
observations  in  the  radius-theta  plane. 
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Fabian  09/02  Number  of  Observations 


Type 

Figure  37.  Number  of  observations  on  2  September  for  each  instrument.  Instrument 
types  are  given  in  Table  2. 
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Figure  38.  Analysis  versus  observations  for  Hurricane  Fabian  on  2  September,  (a) 
Scatterplot  of  all  observations  and  (b)  histogram  of  differences  between 
observations  and  analysis.  Units  of  each  observation  are  the  same  as  in  Table 
5. 
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(a)  (b) 


Figure  39.  Spatial  distribution  of  observations  on  3  September.  Non-radar  observations 
in  the  (a)  radius-height  plane  and  (b)  in  the  radius-theta  plane. 
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Fabian  09/04  Non-radar  observations 


Fabian  09/04  Non-radar  observations 


Figure  40.  Spatial  distribution  of  observations  on  4  September.  Non-radar  observations 
in  the  (a)  radius-height  plane  and  (b)  in  the  radius-theta  plane. 
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(a) 


(b) 


Isabel  09/12  Non-radar  observations 


Figure  41.  Spatial  distribution  of  observations  on  12  September.  Non-radar  observations 
in  the  (a)  radius-height  plane  and  (b)  in  the  radius-theta  plane. 
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Spatial  distribution  of  observations  on  13  September.  Non-radar  observations 
in  the  (a)  radius-height  plane  and  (b)  in  the  radius-theta  plane. 
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Figure  43.  Spatial  distribution  of  observations  on  14  September.  Non-radar  observations 
in  the  (a)  radius-height  plane  and  (b)  in  the  radius-theta  plane. 
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Figure  44.  Fabian  axisymmetric  tangential  wind  (color),  radial  wind  (contour),  and 
secondary  circulation  (vectors)  from  2-4  September. 
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Figure  45. 


Fabian  axisymmetric  angular  momentum  (color),  total  energy  (350  kJ  kg  1  + 
1  kJ  kg’1  contours),  and  secondary  circulation  (vectors)  from  2-4  September. 
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Figure  46.  Isabel  axisymmetric  tangential  wind  (color),  radial  wind  (contour),  and 
secondary  circulation  (vectors)  from  12-14  September. 
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Figure  47.  Isabel  axisymmetric  angular  momentum  (color),  total  energy  (350  kJ  kg  1  +  1 
kJ  kg"1  contours),  and  secondary  circulation  (vectors)  from  12-14  September. 

Table  6.  Statistical  differences  between  analysis  and  all  observations.  Units  of  observations  are 

listed  in  Table  5. 


Mission 

2-Sep 

3-Sep 

4- Sep 

1 2-Sep 

13-Sep 

14-Sep 

#  Obs  xlO5 

7.1105 

8.5040 

5.4717 

9.8731 

6.5533 

17.448 

Minimum 

-42.912 

-37.362 

-29.793 

-53.507 

-67.852 

-43.795 

Maximum 

41.45 

41.632 

36.762 

41.57 

71.787 

34.927 

Mean 

-0.5821 

-0.50257 

-0.29302 

-0.1334 

-0.13643 

-0.3121 

Median 

-0.5448 

-0.4461 

-0.2477 

-0.24437 

-0.069 

-0.35765 

RMS 

5.0665 

4.8648 

5.1959 

5.0483 

6.1022 

4.7359 

Skewness 

0.01482 

0.16971 

0.19735 

0.10709 

0.16345 

0.089939 

Kurtosis 

1.5122 

2.0388 

0.60577 

3.2448 

13.503 

1.4067 

Correlation 

0.99129 

0.99192 

0.98862 

0.99317 

0.99099 

0.99342 
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AIR-SEA  FLUXES 


The  angular  momentum  and  total  energy  budgets  were  resolved  using  the 
SAMURAI-derived  storm  structures  shown  in  the  previous  section.  A  range  of  control 
volume  depths  near  the  top  of  the  inflow  layer  from  400-1000  m  was  utilized,  with  the 
lowest  level  at  both  the  surface  and  at  100  m  height.  Control  volume  widths  from  10-22 
km  in  2  km  increments  were  then  prescribed,  such  that  the  aspect  ratio  of  the  volume 
remained  in  the  20-30  interval  suggested  by  the  sensitivity  analysis.  The  inner  radius  of 
the  volume  was  then  varied  by  2-km  increments  around  a  central  value  of  0.8  normalized 
inner  radius,  which  yields  a  total  of  72  flux  “samples”  per  mission.  An  example  of  the 
control  volumes  used  on  12  September  is  shown  in  Figure  48.  These  volumes  spanned 
the  region  in  and  around  the  eyewall,  such  that  the  average  wind  speed  in  the  volume  was 
always  above  50  m  s’1  for  every  mission.  This  spectrum  of  control  volumes  adequately 
represents  a  reasonable  range  of  shapes  and  sizes  for  obtaining  a  mean  and  standard 
deviation  of  the  budget  residual  at  the  eyewall  on  each  day.  The  exchange  coefficient 
sensitivity  to  the  control  volume  parameters  in  this  range  was  similar  to  that  derived  from 
the  numerical  simulations,  with  the  exception  of  the  WRF  Ck  depth  sensitivity.  A  weak 
decreasing  trend  in  Ck  with  increasing  volume  depth  was  deduced  from  the  observations, 
whereas  no  distinct  decrease  above  400  m  depth  was  found  as  in  the  WRF  simulations. 
Results  with  a  fixed  1  km  control  volume  top  were  within  10%  and  20%  of  the  mean  Cd 
and  Ck  values  over  all  depths,  respectively. 

The  retrieved  surface  stress  values  are  shown  in  Figure  49.  Although  considerable 
spread  exists  in  the  individual  samples,  a  general  agreement  in  the  stress  magnitudes  is 
found  with  a  mean  value  of  9.4  and  standard  deviation  of  4.6  N  m‘  .  A  linear  fit  of  the 
data  has  a  correlation  coefficient  of  0.22  and  only  a  slight  upward  trend  with  wind  speed. 
Missions  on  2  and  13  September  have  the  highest  stress,  but  also  have  two  of  the  deepest 
inflow  layers.  The  14  September  analysis  also  has  a  deep  inflow  layer,  but  with  a  stress 
more  similar  to  the  other  missions.  As  will  be  shown  below,  the  12  September  retrieval 
has  an  anomalously  high  Ck/Cd  ratio  compared  to  the  other  missions,  which  suggests  that 
the  retrieved  stress  is  possibly  too  low.  It  is  not  clear  that  this  retrieval  is  an  outlier,  but 
more  likely  that  the  scatter  is  due  to  inherent  uncertainties  that  were  revealed  in  the 
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sensitivity  analysis.  Removing  this  mission  from  the  sample  improves  the  correlation 
coefficient  to  0.53,  and  results  in  a  more  distinct  upward  trend  with  wind  speed.  Dividing 
by  the  density  integral  and  taking  the  square  root  yields  the  friction  velocity,  which  is 
shown  in  Figure  50.  The  scatter  is  less  than  in  the  stress  retrieval  due  to  the  square  root 
dependence,  with  a  mean  value  of  2.8  m  s"1  and  a  0.7  m  s'1  standard  deviation.  These 
estimates  also  indicate  a  slight  upward  trend  with  wind  speed  with  a  correlation  of  0.24, 
or  0.53  with  the  12  September  mission  removed. 

The  retrieved  enthalpy  flux  is  shown  in  Figure  51.  These  are  the  first  known 
estimates  of  the  enthalpy  flux  at  wind  speeds  greater  than  50  m  s’1,  and  likely  the  highest 
estimates  of  heat  flux  recorded  over  the  ocean  surface.  There  is  less  scatter  than  in  the 
stress  retrieval,  and  an  increasing  linear  correlation  with  wind  speed  of  0.81.  While  some 
of  the  individual  samples  clearly  are  too  low  (near  zero  on  3  September)  and  too  high 
(4000  W  m’  on  12  September),  the  mean  enthalpy  flux  for  each  mission  and  trend 
appears  reasonable,  with  a  variation  from  764  W  m’~  at  52  m  s'  to  2189  W  m"  at 
72  m  s'1.  The  reasons  for  the  reduced  scatter  and  better  correlation  in  the  energy  budget 
compared  to  the  momentum  budget  are  not  apparent  but,  given  the  high  variability  in  the 
fluxes  from  numerical  model  simulations,  this  result  also  is  not  too  disconcerting.  It 
would  appear  that  the  derived  thermodynamic  structure  was  steadier  on  the  ~6  hour 
mission  timescale  than  the  kinematic  structure,  but  it  is  impossible  to  validate  this  with 
the  current  dataset. 

The  derived  drag  and  enthalpy  exchange  coefficients  are  shown  in  Figure  52.  The 
relatively  low  scatter  of  both  results  is  very  good  considering  the  potential  errors,  and  the 
enthalpy  exchange  coefficients  have  a  better  agreement  than  do  the  drag  coefficients,  as 
would  be  expected  from  the  previous  retrievals.  The  mean  drag  coefficient  is  2.4  xlO' 
with  a  standard  deviation  of  1.1  xlO'  ,  which  corresponds  to  ±16%  uncertainty  at  the  67% 
confidence  level,  and  a  d9 3%  uncertainty  at  the  95%  confidence  level.  The  mean 
enthalpy  exchange  coefficient  is  1.0  xlO"^  with  a  standard  deviation  of  0.4  xlO"  ,  which  is 
a  ±10%  uncertainty  at  the  67%  confidence  level,  and  an  d80%  uncertainty  at  the  95% 
confidence  level.  These  percentage  uncertainties  are  consistent  with  the  sensitivity 
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analysis  in  Chapter  III  that  indicated  a  -50%  standard  deviation  primarily  associated  with 
the  neglected  budget  terms.  The  drag  coefficient  results  do  not  indicate  a  significant  wind 
speed  dependence,  with  only  a  slightly  decreasing  linear  correlation  of  only  0.11.  The 
enthalpy  coefficient  shows  a  slight  increase  with  wind  speed,  with  a  slope  of 
0.03  (m  s"1)'1  and  a  higher  linear  correlation  of  0.53.  Given  the  uncertainties  and  the  weak 
slopes  of  both  fits,  there  is  no  statistically  significant  change  in  either  the  drag  or  enthalpy 
coefficient  in  this  wind  speed  range. 

Histograms  of  the  Cd  and  Ck  samples  (Figure  53)  suggest  that  the  distribution  of 
the  samples  is  similar  and  approximately  normal  in  both  cases.  Low  right  skewness  (0.56 
and  0.45  for  Cd  and  Ck,  respectively)  with  a  slightly  rounder  peak  than  normal  and 
negative  kurtosis  (-0.58  and  -0.12,  respectively)  are  found  for  both  coefficients.  These 
distributions  yield  good  confidence  that  the  sample  mean  values  of  the  distribution  are 
representative  of  the  actual  mean  bulk  exchange  coefficients  at  these  wind  speeds,  with 
the  majority  of  the  variability  attributable  to  random  errors  associated  with  the 
uncertainties  examined  in  Chapter  III.  The  positive  skewness  suggests  that  there  may  be  a 
small  low  bias  in  the  mean,  but  this  uncertainty  is  well  within  the  standard  deviation  of 
the  samples.  The  cumulative  probability  distribution  (Figure  53,  right)  graphically 
illustrates  the  probability  that  the  exchange  coefficients  are  at  or  below  a  particular 
magnitude.  The  drag  coefficient  has  a  steep  slope  above  1  xlO"3,  with  -50%,  70%,  and 
90%  probabilities  that  Cd  is  less  than  2,  3,  and  4  xlO"  ,  respectively.  The  maximum 
estimated  drag  coefficient  is  less  than  5  xlO^  with  99%  probability,  although  this  outlier 
estimate  is  not  that  useful  given  prior  estimates’  magnitude  (Chapter  II.B).  The  slope  of 
the  Ck  distribution  is  much  flatter  than  the  drag  coefficient  distribution.  The  maximum 
estimated  enthalpy  exchange  coefficient  is  2  xlO"  with  99%  probability,  with  an  -60% 
probability  that  the  value  is  less  than  1  xlO'  .  These  probabilities  assume  that  the  432 
samples  obtained  from  varying  the  control  volume  over  each  of  the  six  missions 
accurately  represent  samples  from  the  true  population  distribution  of  exchange  coefficient 
magnitudes.  Since  the  numerical  simulation  sensitivity  tests  indicated  that  the 
uncertainties  associated  with  these  control  volumes  may  be  similar  to  those  obtained  with 


112 


the  real  data,  and  the  other  error  sources  are  also  primarily  random  over  the  six  missions, 
this  is  believed  to  be  a  fair  assumption. 

Dividing  the  mean  Ck/Cd  gives  an  average  ratio  of  0.4,  with  a  variation  from  0.17 
to  1.05  by  adding  and  subtracting  a  single  standard  deviation  to  each  coefficient  in  each 
direction.  The  ratios  of  Ck/Cd  calculated  from  the  different  samples  and  the  cumulative 
probability  distribution  are  shown  in  Figure  54.  With  the  exception  of  12  September,  the 
retrieved  ratios  are  in  relatively  good  agreement.  The  mean  value  of  all  individually 
calculated  ratios  is  0.48  with  a  standard  deviation  of  0.27,  which  is  similar  to  the  0.4  ratio 
and  uncertainty  obtained  by  dividing  the  mean  Ck  and  Cd  from  all  the  samples.  A  linear 
fit  shows  an  increasing  ratio  with  wind  speed  with  a  correlation  of  0.58,  but  this  is 
primarily  due  to  the  large  ratios  above  70  m  s'1.  It  is  unclear  why  the  12  September 
mission  appears  to  be  an  outlier,  but  this  is  probably  due  to  unresolved  errors  as  opposed 
to  a  significant  change  in  the  ratio  above  70  m  s'1.  It  cannot  be  ruled  out  that  the  drag 
coefficient  is  reduced  at  these  wind  speeds.  Since  the  enthalpy  exchange  coefficient  for 
this  mission  agrees  well  with  the  other  five,  the  drag  coefficient  derived  for  this  day 
probably  is  too  low,  although  this  discrepancy  could  be  because  Hurricane  Isabel  may 
have  undergone  an  eyewall  replacement  cycle  shortly  after  this  mission,  and  this  cycle 
was  accompanied  with  a  change  in  the  angular  momentum  structure  and  a  larger 
integrated  budget  tendency  term.  Even  with  this  mission  included,  the  cumulative 
probability  distribution  of  all  Ck/Cd  ratios  indicates  an  -80%  probability  that  the  ratio  is 
less  than  0.75,  and  a  93%  probability  it  is  less  than  1.0.  Without  the  12  September  values 
these  probabilities  increase  to  91%  and  97%,  respectively  (not  shown). 
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Figure  48. 


Control  volumes  used  for  flux  retrieval  on  12  September.  Tangential  wind 
(color)  and  secondary  circulation  (vectors)  are  shown  for  reference. 


Surface  Wind  Speed  (m  s  ') 


Figure  49.  Derived  surface  stress  (T-0)  from  budget  retrieval.  Gray  dots  indicate 
individual  samples  from  different  control  volumes,  and  large  symbols 
indicate  mean  values  from  each  research  mission.  Error  bars  indicate  one 
standard  deviation  in  Tze  (vertical)  and  average  surface  wind  speed 
(horizontal). 
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Friction  Velocity 
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Figure  50.  As  in  Figure  49,  except  derived  friction  velocity  (u*)  from  budget  retrievals. 
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As  in  Figure  49,  except  derived  enthalpy  flux  from  budget  retrievals 
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Figure  5 1 . 


(a)  Drag  Coefficient 
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(b)  Enthalpy  Exchange  Coefficient 
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Figure  52.  As  in  Figure  49,  except  derived  bulk  exchange  coefficients  from  budget 
retrievals. 
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Figure  53.  Flistogram  of  the  percentage  of  samples  in  0.5xl0'3  bins  for  Cd  (black)  and 
Ck  (gray)  for  all  retrievals  (left).  Cumulative  probability  distributions  for  Cd 
(black)  and  Ck  (gray)  retrievals  (right). 
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Figure  54.  Ratio  of  Ck/Cd  from  budget  method.  Left  panel  has  symbols  as  in  Figure  49. 

Right  panel  shows  cumulative  probability  distribution  from  all  samples. 
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VI.  SUMMARY  AND  CONCLUSIONS 


The  magnitude  of  surface  fluxes  in  major  hurricanes  has  been  very  difficult  to 
determine,  and  the  resulting  uncertainty  remains  an  important  aspect  of  both  research  and 
operations  involving  tropical  cyclones.  For  the  first  time,  enthalpy  fluxes,  the  bulk 
enthalpy  exchange  coefficient,  and  the  ratio  of  Ck/Cd  have  been  estimated  at  major 
hurricane  wind  speeds  using  a  budget  methodology.  The  results  presented  in  this  study 
also  augment  the  limited  stress  and  drag  coefficient  estimates  above  50  m  s'1.  This  study 
improves  upon  previous  studies  using  a  budget  method  for  estimating  the  surface  fluxes 
in  hurricanes  by  expanding  the  budgets  to  include  total  energy,  by  using  improved  data 
quality  and  quantity,  by  the  use  of  a  more  sophisticated  analysis  technique,  and  by 
including  a  comprehensive  error  analysis.  The  six  CBLAST  missions  into  major 
Hurricanes  Fabian  and  Isabel  provided  a  unique  opportunity  to  apply  this  budget  method 
at  a  range  of  surface  wind  speeds  from  52-72  m  s'1  and  to  obtain  a  good  statistical 
sampling  of  flux  retrievals.  The  new  variational  analysis  scheme  called  SAMURAI  was 
developed  from  first  principles  to  determine  the  maximum  likelihood  estimate  of  the 
axisymmetric  structure  for  the  budgets  by  incorporating  flight  level,  dropsonde,  SFMR, 
and  Doppler  radar  data.  Although  the  derived  exchange  coefficients  have  some 
uncertainty,  the  sensitivity  analysis  using  simulated  data  has  allowed  for  a  mitigation  of 
some  of  the  very  large  potential  errors,  and  provides  confidence  in  the  quantitative 
uncertainty  estimates.  The  main  sources  of  error  identified  were:  residuals  in  the  mass 
continuity  equation;  unresolved  volume-integrated  tendency  terms;  unresolved  vertical 
eddy  flux  terms  at  the  top  of  control  volume;  sea- surface  temperature  errors;  circulation 
center  errors;  and  analysis  errors.  The  “shear  flux”  term  also  was  found  to  be  important  in 
the  energy  budget,  and  was  included  through  an  estimate  of  the  drag  coefficient  obtained 
from  the  angular  momentum  budget.  A  total  of  432  samples  from  72  control  volumes 
over  the  six  analysis  days  provides  a  good  sampling  of  the  data  that  takes  into  account 
errors  associated  with  the  budget  retrieval. 

Estimates  of  the  drag  coefficient  from  the  current  research  in  relation  to  previous 
studies  are  shown  in  Figure  55.  The  black  symbols  are  taken  from  French  et  al.  (2007) 
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that  represent  the  eddy  correlation  CBLAST  estimates  and  the  laboratory  estimates  from 
Donelan  et  al.  (2004).  The  blue  symbols  are  taken  from  Vickery  et  al.  (2009)  showing  the 
flux-profile  estimates  obtained  from  dropsondes.  The  current  results  are  shown  in  green, 
along  with  95%  confidence  intervals.  Although  the  spread  in  the  budget  estimates  from 
this  study  is  non-trivial,  the  retrieved  magnitudes  are  in  general  agreement  with  the 
previous  studies  that  the  drag  coefficient  does  not  continue  to  increase  beyond  ~30  ms"1. 
At  the  upper  end  of  the  95%  confidence  level,  the  mean  drag  coefficient  from  some  of  the 
runs  is  near  the  extrapolated  Large  and  Pond  (1981)  curve  (dash-dot  black  curve),  but  the 
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mean  drag  coefficient  averaged  over  all  the  samples  is  lower  at  2.4  xlO"  with  a  ±16% 
uncertainty  at  the  67%  confidence  level.  At  the  edge  of  the  probability  distribution,  the 
93%  uncertainty  at  the  95%  confidence  level  does  not  preclude  the  possibility  that  the 
drag  coefficient  is  near  the  extrapolated  value,  but  the  cumulative  distribution  indicates 

■j 

that  the  magnitude  has  a  90%  probability  of  being  less  than  4  xl0"\ 

The  individual  enthalpy  exchange  coefficient  samples  are  shown  with  previous 
estimates  in  Figure  56.  A  similar  degree  of  scatter  exists  in  the  present  results  as  those 
from  the  HEXOS  and  CBLAST  eddy  correlations,  and  the  mean  value  of  1.0  xlO"  with 
an  d80%  uncertainty  at  the  95%  confidence  level  is  very  similar  to  both  the  field  and 
laboratory  results.  The  good  0.81  linear  correlation  of  the  enthalpy  flux  with  wind  speed, 
and  consistency  of  the  exchange  coefficients  from  the  different  missions  yields  good 
confidence  in  these  results.  Viewed  as  percentage  errors,  the  uncertainty  in  the  drag  and 
the  enthalpy  exchange  coefficients  is  similar.  However,  the  enthalpy  coefficient  has  a 
lower  uncertainty  in  terms  of  absolute  error.  A  linear  fit  of  the  enthalpy  coefficient 
indicates  a  slight  increase  with  wind  speed  above  50  m  s"1,  but  this  increase  is  well  within 
the  uncertainty  range  and  cannot  be  concluded  definitively.  This  statistical  sample 
indicates  that  the  enthalpy  coefficient  has  an  -70%  probability  of  being  equal  to  or  less 
than  1.2x10"  ,  which  is  the  approximate  value  determined  at  wind  speeds  greater  than 
15  ms"1.  These  results  suggest  that  it  is  probable  that  the  enthalpy  exchange  coefficient  is 
not  dependent  on  wind  speed  in  hurricane  conditions.  Since  the  current  results  implicitly 
include  the  effects  of  sea  spray,  this  would  also  suggest  that  spray  effects  do  not  change 
the  enthalpy  exchange  coefficient.  One  hypothesis  is  that  the  spray  flux  has  simply 
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replaced  the  interfacial  flux  at  these  wind  speeds,  which  leads  to  a  similar  overall 
enthalpy  flux.  However,  this  cannot  be  assessed  by  the  current  methodology. 

The  results  from  this  study  have  several  theoretical  and  practical  implications. 
The  magnitudes  of  the  drag  and  enthalpy  exchange  coefficients  above  50  m  s'1  provide 
some  support  of  the  similarity  hypothesis  of  Emanuel  (2003)  that  the  drag  and  enthalpy 
coefficients  would  remain  constant  at  extreme  wind  speeds.  The  retrieved  Cr/Cd  ratios 
from  the  present  study  are  shown  with  previous  estimates  in  Figure  57.  These  results 
indicate  that  the  ratio  is  likely  less  than  the  -1.0  estimate  derived  by  Emanuel  (2003)  at 
the  -28  C  temperature  range,  and  perhaps  may  be  as  low  as  0.4.  The  uncertainty  in  the 
coefficients  and  limited  sea-surface  temperature  range  precludes  testing  of  the  hypothesis 
that  the  enthalpy  exchange  coefficient  is  SST-dependent.  There  is  an  -80%  probability 
from  these  retrievals  that  the  ratio  is  less  than  the  0.75  ratio  proposed  as  a  threshold  for 
hurricane  development  by  Emanuel  (1995).  Recent  theoretical  and  computational 
analyses  suggest  that  this  discrepancy  is  largely  due  to  gradient  wind  imbalance  in  the 
boundary  layer  (Smith  et  al.  2008;  Bryan  and  Rotunno  2009).  These  results  are  also 
consistent  with  recent  three-dimensional  numerical  simulations  that  intensified  to  major 
hurricane  status  with  Cr/Cd  ratios  as  low  as  0.1  (Montgomery  et  al.  2010).  Although 
evidence  is  accumulating  that  potential  intensity  may  not  be  as  sensitive  to  Cr/Cd  as 
originally  formulated,  this  does  not  diminish  the  importance  of  having  accurate  estimates 
of  these  parameters  at  major  hurricane  wind  speeds.  These  results  suggest  that  the  lower 
bound  of  the  E-PI  range  for  Hurricane  Isabel  presented  in  BM08  would  be  the  most 
accurate  (c.f.,  their  Figure  1),  which  yields  potential  intensity  estimates  below  50  m  s'1. 
These  new  estimates  of  Cr  and  Cd  should  help  to  improve  potential  intensity  theory  and 
understanding  of  tropical  cyclone  intensity  change.  From  a  modeling  perspective,  these 
estimates  should  provide  additional  confidence  in  an  improved  physical  basis  for  surface 
layer  schemes  that  do  not  increase  the  drag  or  enthalpy  coefficient  as  the  wind  speeds 
increase  above  50  m  s'1.  A  continued  examination  of  the  physics  parameterizations  at 
high  wind  speeds  in  WRF  and  COAMPS,  especially  the  new  COAMPS-TC  that  is 
currently  under  development,  will  be  conducted  as  part  of  future  research. 
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Several  additional  areas  beyond  this  study  deserve  further  attention.  With  this 
proof  of  concept  in  place,  fluxes  in  other  storms  with  sufficient  research  data  quality 
could  be  examined  to  improve  the  statistical  sampling  and  potentially  expanding  the  wind 
speed  ranges  presented  in  this  paper.  Several  tropical  cyclones  with  good  reconnaissance 
and  research  data  exist  that  would  be  candidates.  In  order  to  improve  the  method,  the 
budget  error  analysis  indicates  that  two  of  the  largest  sources  of  error  are  the  unresolved 
tendency  term  and  circulation  center  errors.  One  possible  approach  would  be  to  relax  the 
axisymmetric  constraint,  and  reformulate  the  budgets  in  a  Cartesian  coordinate  system 
centered  on  a  particular  region  of  the  storm.  This  would  remove  truncation  errors 
associated  with  the  axisymmetric  cylindrical  transform,  but  would  likely  increase  errors 
associated  with  the  unresolved  terms.  If  multiple  passes  could  be  conducted  through  a 
region,  the  tendency  term  could  be  estimated  on  relatively  short  timescales.  The  inclusion 
of  radar  data  to  completely  document  the  kinematics  in  the  control  volume  would  be 
optimal.  Since  it  is  unclear  whether  this  approach  would  reduce  the  uncertainty,  a 
comprehensive  sensitivity  analysis  similar  to  the  one  conducted  herein  is  recommended. 
Lastly,  the  SAMURAI  analysis  package  will  continue  to  be  developed,  with  a  natural 
extension  to  three  dimensions  and  the  Cartesian  domain,  and  with  the  addition  of  other 
data  sources  and  analysis  variables  that  are  relevant  for  mesoscale  studies.  The  usefulness 
of  the  software  for  tropical  cyclone  data  assimilation  applications  with  numerical  weather 
models  also  will  be  investigated. 


Figure  55.  Mean  drag  coefficients  from  this  study  (green  circles)  compared  with 
previous  studies.  Black  symbols  adapted  from  French  et  al.  (2007)  and  blue 
symbols  adapted  from  Vickery  et  al.  (2009).  Red  line  indicates  measured 
(thick)  and  extrapolated  (thin)  Large  and  Pond  (1981)  drag  coefficient. 
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Figure  56.  Wind  speed  dependence  of  Ck  from  this  study  (green  squares)  compared  with 
previous  studies  as  summarized  by  Haus  et  al.  (2010).  ASIST  laboratory 
results  (blue  circles)  and  CBLAST  (red  triangles)  measurements  shown  with 
HEXOS  results  (gray  x’s).  The  mean  and  95%  confidence  intervals  are 
shown  in  black. 
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Figure  57.  Wind  speed  dependence  of  Ck/Cd  from  this  study  (green  squares)  compared 
with  previous  studies  as  summarized  by  Flaus  et  al.  (2010).  ASIST  laboratory 
results  (blue  circles)  and  CBLAST  (red  triangles)  measurements  shown  with 
HEXOS  results  (gray  x’s).  The  mean  and  95%  confidence  intervals  are 
shown  in  black.  Purple  dashed  line  indicates  0.75  ratio. 


123 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


124 


LIST  OF  REFERENCES 


Aberson,  S.  D.,  M.  T.  Montgomery,  M.  Bell,  and  M.  Black,  2006:  Superintense  winds  in 
hurricane  Isabel  (2003).  Part  II:  Extreme  wind  speeds.  Bull.  Amer.  Met.  Soc.,  87, 
1349-1354. 

Atlas,  D.,  R.  C.  Srivastava,  and  R.  S.  Sekhon,  1973:  Doppler  radar  characteristics  of 
precipitation  at  vertical  incidence.  Rev.  Geophys.  Space  Phys.,  11,  1-35. 

Bao,  J.-W.,  S.  A.  Michelson,  and  J.  M.  Wilczak,  2002:  Sensitivity  of  numerical 

simulations  to  parameterizations  of  roughness  for  surface  heat  fluxes  at  high 
winds  over  the  sea.  Mon.  Wea.  Rev.,  130,  1926-1932. 

Bargen,  D.  W.,  and  R.  C.  Brown,  1980:  Interactive  radar  velocity  unfolding.  Preprints, 
19th  Conference  on  Radar  Meteorology,  Miami  Beach,  Amer.  Meteor.  Soc.,  278- 
285. 

Barnes,  S.  L.,  1973:  Mesoscale  objective  analysis  using  weighted  time-series 

observations.  NOAA  Tech.  Memo.  ERL  NSSL-62,  National  Severe  Storms 
Laboratory,  Norman,  OK  73069,  60  pp. 

Beard,  K.,  1985:  Simple  altitude  adjustments  to  raindrop  velocities  for  Doppler  radar 
analysis,  J.  Atmos,  and  Oceanic  Technol.,  2,  468-471 

Bell,  M.  M.,  2006:  Observed  structure,  evolution,  and  potential  intensity  of  category  five 
Hurricane  Isabel  (2003)  from  12-14  September.  M.S.  thesis,  Dept,  of 
Atmospheric  Science,  Colorado  State  University,  98  pp. 

- ,  and  M.  T.  Montgomery,  2008:  Observed  structure,  evolution,  and  potential 

intensity  of  category  five  Hurricane  Isabel  (2003)  from  12-14  September.  Mon. 
Wea.  Rev.,  136,  2023-2045. 

Black,  M.  L.,  R.  W.  Burpee,  and  F.  D.  Marks  Jr.,  1996:  Vertical  motion  characteristics  of 
tropical  cyclones  determined  with  airborne  Doppler  radial  velocities.  J.  Atmos. 
Sci.,  53,  1887-1909. 

Black,  M.  L.,  E.  A.  D’Asaro,  W.  M.  Drennan,  J.  R.  French,  P.  P.  Niiler,  T.  B.  Sanford,  E. 
J.  Terrill,  E.  J.  Walsh,  and  J.  Zhan,  2007:  Air-sea  exchange  in  hurricanes: 
Synthesis  of  observations  from  the  Coupled  Boundary  Layer  Air-Sea  Transfer 
experiment.  Bull.  Amer.  Met.  Soc.,  88,  357-374. 

Bosart,  B.  L.,  W.-C.  Lee,  and  R.  M.  Wakimoto,  2002:  Procedures  to  improve  the 

accuracy  of  airborne  Doppler  radar  data.  J.  Atmos.  Oceanic  Technol.,  19,  322- 
339. 


125 


Braun,  S.  A.,  and  W.-K.  Tao,  2000:  Sensitivity  of  high-resolution  simulations  of 

hurricane  Bob  (1991)  to  planetary  boundary  layer  parameterizations.  Mon.  Wea. 
Rev.,  128,  3941-3961. 

Bryan,  G.  H.,  and  R.  Rotunno,  2009:  The  maximum  intensity  of  tropical  cyclones  in 
axisymmetric  numerical  model  simulations.  Mon.  Wea.  Rev.,  137,  1770-1789. 

Courtier,  P.,  J.-N.  Thepaut,  and  A.  Hollingsworth,  1994:  A  strategy  for  operational 
implementation  of  4D-VAR,  using  an  incremental  approach.  Quart.  J.  Roy. 
Meteor.  Soc,  120,  1367-1387. 

Cram,  T.  A.,  J.  Persing,  M.  T.  Montgomery,  and  S.  A.  Braun,  2005:  A  Lagrangian 

trajectory  view  on  transport  and  mixing  processes  between  the  eye,  eyewall  and 
environment  using  a  high  resolution  simulation  of  Hurricane  Bonnie  (1998).  J. 
Atmos.  Sci.,  64,  1835-1856. 

Cressman,  G.  P.,  1959:  An  operational  objective  analysis  system.  Mon.  Wea.  Rev.,  87, 
367-374. 

Davis,  C.,  W.  Wang,  S.  Chen,  Y.  Chen,  K.  Corbosiero,  M.  DeMaria,  J.  Dudhia,  G. 
Holland,  J.  Klemp,  J.  Michalakes,  H.  Reeves,  R.  Rotunno,  C.  Snyder,  and  Q. 

Xiao,  2008:  Prediction  of  landfalling  hurricanes  with  the  Advanced  Hurricane 
WRF  model.  Mon.  Wea.  Rev.,  136,  1990-2005. 

Deardorff,  J.  W.,  1972:  Parameterization  of  the  planetary  boundary  layer  for  use  in 
general  circulation  models.  Mon.  Wea.  Rev.,  100,  93-106. 

DeCosmo,  J.,  K.  B.  Katsaros,  S.  D.  Smith,  R.  J.  Anderson,  W.  A.  Oost,  K.  Bumke,  and  H. 
Chadwick,  1996:  Air-sea  exchange  of  water  vapor  and  sensible  heat:  The 
Humidity  Exchange  over  the  Sea  (HEXOS)  results. ./.  Geophys.  Res.,  101,  12001- 
12016. 

Donelan,  M.  A.,  B.  K.  Haus,  N.  Reul,  W.  J.  Plant,  M.  Stianssnie,  H.  C.  Graber,  O.  B. 

Brown,  and  E.  S.  Saltzman,  2004:  On  the  limiting  aerodynamic  roughness  of  the 
ocean  in  very  strong  winds.  Geophys.  Res.  Lett.,  31,  LI 8306. 

Drennan,  W.  M.,  J.  A.  Zhang,  J.  R.  French,  C.  McCormick,  and  P.  G.  Black,  2007: 
Turbulent  fluxes  in  the  hurricane  boundary  layer.  Part  II:  Latent  heat  flux.  J. 
Atmos.  Sci.,  64,  1103-1115. 

Eastin,  M.  D.,  P.  G.  Black,  and  W.  M.  Gray,  2002:  Flight-level  instrument  wetting  errors 
in  hurricanes.  Part  I:  Observations.  Mon.  Wea.  Rev.,  130,  825-841.380. 

Emanuel,  K.A.,  1986:  An  air-sea  interaction  theory  for  tropical  cyclones.  Part  I:  Steady- 
state  maintenance.  J.  Atmos.  Sci.  43,  585-604. 


126 


- ,  1995:  Sensitivity  of  tropical  cyclones  to  surface  exchange  coefficients  and  a 

revised  steady-state  model  incorporating  eye  dynamics.  J.  Atmos.  Sci.,  52,  3969- 
3976. 

- ,  2003:  A  similarity  hypothesis  for  air-sea  exchange  at  extreme  wind  speeds.  J. 

Atmos.  Sci.,  60,  1420-1428. 

French,  J.  R.,  W.  M.  Drennan,  J.  A.  Zhang,  and  P.  G.  Black,  2007:  Turbulent  fluxes  in  the 
hurricane  boundary  layer.  Part  I:  Momentum  flux.  J.  Atmos.  Sci.,  64,  1089-1102. 

Gao,  J.,  M.  Xue,  K.  Brewster,  and  K.  K.  Droegemeier,  2004:  A  three-dimensional 

variational  data  analysis  method  with  recursive  filter  for  Doppler  radars.  J.  Atmos. 
Oceanic  Technol,  21,  457-469. 

Haus,  B.  K.,  D.  Jeong,  M.  A.  Donelan,  J.  A.  Zhang,  and  I.  Savelyev,  2010:  Relative  rates 
of  sea-air  heat  transfer  and  frictional  drag  in  very  high  winds.  Geophys.  Res.  Lett., 
37,  L07802,  doi:10.1029/2009GL042206. 

Hausman,  S.  A.,  K.  V.  Ooyama,  and  W.  H.  Schubert,  2006:  Potential  vorticity  structure 
of  simulated  hurricanes.  J.  Atmos.  Sci.,  63,  87-108. 

Hawkins,  H.  F.,  and  D.  T.  Rubsam,  1968:  Hurricane  Hilda,  1964.  II.  Structure  and 
budgets  of  the  hurricane  on  October  1,  1964.  Mon.  Wea.  Rev.  96,  617-636. 

- ,  and  S.  M.  Imbembo,  1976:  The  structure  of  a  small,  intense  hurricane — Inez  1966. 

Mon.  Wea.  Rev.  104,  418-442. 

Hayden,  C.  M.,  and  R.  J.  Purser,  1995:  Recursive  filter  for  objective  analysis  of 

meteorological  fields:  Applications  to  NESDIS  operational  processing.  J.  Appl. 
Meteor,  34,  3-15. 

Hong,  S.-Y.,  Y.  Noh,  and  J.  Dudhia,  2006:  A  new  vertical  diffusion  package  with  an 
explicit  treatment  of  entrainment  processes.  Mon.  Wea.  Rev.,  134,  2318-2341. 

Huang,  X.-Y.,  2000:  Variational  analysis  using  spatial  filters.  Mon.  Wea.  Rev,  128,  2588- 
2600. 

Hock,  T.  F.,  and  J.  L.  Franklin,  1999:  The  NCAR  GPS  dropwindsonde.  Bull.  Amer. 
Meteor.  Soc.,  80,  407-420. 

Hughes,  L.  A.,  1952:  On  the  low-level  wind  structure  of  tropical  storms.  J.  Meteor.,  9, 
422-428. 

Jarosz,  E.,  D.  A.  Mitchell,  D.  W.  Wang,  and  W.  J.  Teague,  2007:  Bottom-up 

determination  of  air-sea  momentum  exchange  under  a  major  tropical  cyclone. 
Science,  315,  1707-1709. 


127 


Jeong,  D.,  B.  K.  Haus,  and  M.  E.  Donelan,  2010:  Enthalpy  transfer  across  the  air-water 
interface  in  high  winds  including  spray.  J.  Atmos.  Sci.,  submitted. 

Jordan,  C.  L.,  1952:  On  the  low-level  structure  of  the  typhoon  eye.  J.  Meteor.,  9,  285- 
290. 

- ,  1958:  Mean  soundings  for  the  West  Indies  area.  J.  Meteor.,  15,  91-97. 

Joss,  J.,  and  A.  Waldvogel,  1970:  A  method  to  improve  the  accuracy  of  radar-measured 
amounts  of  precipitation.  Proc.  14th  Radar  Meteorology  Conf.,  Tucson,  AZ, 
Amer.  Meteor.  Soc.,  237-238. 

Keeler,  R.  J.,  and  S.  M.  Ellis,  2000:  Observational  error  covariance  matrices  for  radar 
data  assimilation.  Phys.  Chem.  Earth,  25,  1277-1280. 

Kessler,  E.,  1969:  On  distribution  and  continuity  of  water  substance  in  atmospheric 
circulations.  American  Meteorological  Society,  Meteorol.  Monogr.,  10,  32-84. 

Kleinschmidt,  E.,  1951:  Grundlagen  einer  theorie  der  tropischen  zyklonen.  Arch. 
Meteorol.  Geophys.  Bioklimatol.,  A4,  53-72. 

Large,  W.  G.,  and  S.  Pond,  1981:  Open  ocean  momentum  flux  measurements  in  moderate 
to  strong  winds.  J.  Phys.  Oceanogr.,  11,  324-336. 

- ,  and - ,,  1982:  Sensible  and  latent  heat  flux  measurements  over  the  ocean.  J. 

Phys.  Oceanogr.,  12,  464-482. 

Lee,  W-C.,  and  F.  D.  Marks,  2000:  Tropical  cyclone  kinematic  structure  retrieved  from 
single  Doppler  radar  observations.  Part  II:  The  GBVTD-simplex  center  finding 
algorithm.  Mon.  Wea.  Rev.,  128,  1925-1936. 

Lorenc,  A.  C.,  1986:  Analysis  methods  for  numerical  weather  prediction.  Quart.  J.  Roy. 
Meteor.  Soc,  112,  1 177-1 194. 

Malkus,  J.  S.,  and  H.  Riehl,  1960:  On  the  dynamics  and  energy  transformations  in  steady- 
state  hurricanes.  Tellus,  12,  1-20. 

Marks,  F.  D.  Jr.,  R.  A.  Houze,  and  J.  Gamache,  1992:  Dual-aircraft  investigation  of  the 
inner  core  of  Hurricane  Norbert.  Part  I:  Kinematic  structure.  J.  Atmos.  Sci.,  49, 
919-942. 

Miller,  B.  I.,  1962:  “On  the  momentum  and  energy  balance  of  Hurricane  Helene  (1958),” 
National  Hurricane  Research  Project  Report  No.  53,  U.S.  Weather  Bureau, 
Washington,  D.C.,  19  pp. 

- ,  1964:  A  study  of  the  filling  of  Hurricane  Donna  1960:  Over  land.  Mon.  Wea.  Rev., 

92,  389-406. 


128 


- ,  1969:  Experiments  in  forecasting  hurricane  development  with  real  data.  National 

Hurricane  Research  Laboratory  Report,  Miami,  Fla.,  28  pp. 

Mohr,  C.  G.,  1988:  CEDRIC — Cartesian  Space  Data  Processor.  National  Center  for 
Atmospheric  Research,  Boulder,  CO,  78  pp. 

Montgomery,  M.  T.,  M.  M.  Bell,  S.  D.  Aberson,  and  M.  L.  Black,  2006:  Hurricane  Isabel 
(2003):  New  insights  into  the  physics  of  intense  storms.  Part  I:  Mean  vortex 
structure  and  maximum  intensity  estimates.  Bull.  Amer.  Meteor.  Soc .,  87,  1335- 
1347. 

- ,  R.  Smith,  and  S.  Nguyen,  2010:  Sensitivity  of  tropical  cyclone  models  to  the 

surface  drag  coefficient.  Quart.  J.  Roy.  Meteor.  Soc,  submitted. 

Moss,  M.  S.,  and  S.  L.  Rosenthal,  1975:  On  the  estimation  of  planetary  boundary  layer 
variables  in  mature  hurricanes.  Mon.  Wea.  Rev ,  103,  980-988. 

Nelder,  J.  A.,  and  R.  Mead,  1965:  A  simplex  method  for  function  minimization.  Comput. 
J.,7,  308-313. 


Noh,  Y.,  W.  G.  Cheon,  S.  Y.  Hong,  and  S.  Raasch,  2003:  Improvement  of  the  K-profile 
model  for  the  planetary  boundary  layer  based  on  large  eddy  simulation  data. 
Bound. -Layer  Meteor.,  107,  421-427. 

Nolan,  D.  S.,  M.  T.  Montgomery,  and  L.  D.  Grasso,  2001:  The  wavenumber-one 

instability  and  trochoidal  motion  of  hurricane-like  vortices.  J.  Atmos.  Sci,  58, 
3243-3270. 

Ooyama,  K.  V.,  1969:  Numerical  simulation  of  the  life  cycle  of  tropical  cyclones.  J. 
Atmos.  Sci.,  26,  3-40. 

- ,  1987:  Scale  controlled  objective  analysis.  Mon.  Wea.  Rev.,  115,  2479-2506. 

- ,  2002:  The  cubic-spline  transform  method:  Basic  definitions  and  tests  in  a  ID  single 

domain.  Mon.  Wea.  Rev.,  130,  2392-2415. 

Palmen,  E.,  and  H.  Riehl,  1957:  Budget  of  angular  momentum  and  energy  in  tropical 
cyclones.  J.  Meteor,  14,  150-159. 

Persing,  J.,  and  M.  T.  Montgomery,  2003:  Hurricane  superintensity.  J.  Atmos.  Sci.  60, 
2349-2371. 

Powell,  M.  D.,  P.  J.  Vickery,  and  T.  A.  Reinhold,  2003:  Reduced  drag  coefficient  for 
high  wind  speeds  in  tropical  cyclones.  Nature,  422,  279-283. 

Press,  W.  H.,  S.  A.  Teukolsky,  W.  T.  Vetterling,  and  B.  P.  Flannery,  2002:  Numerical 

recipes  in  C++,  the  art  of  scientific  computing.  Cambridge  Univ.  Press,  1002  pp. 


129 


Purser,  R.  J.,  W.-S.  Wu,  D.  Parrish,  and  N.  M.  Roberts,  2003:  Numerical  aspects  of  the 
application  of  recursive  filters  to  variational  statistical  analysis.  Part  I:  Spatially 
homogeneous  and  isotropic  Gaussian  covariances.  Mon.  Wea.  Rev,  131,  1524— 
1535. 

Rappaport,  E.  N.,  J.  Franklin,  L.  Avila,  S.  Baig,  J.  Beven,  E.  Blake,  C.  Burr,  J.-G.  Jiing, 
C.  Juckins,  R.  Kanbb,  C.  Landsea,  M.  Mainelli,  M.  Mayfield,  C.  McAdie,  R. 
Pasch,  C.  Cisko,  S.  Stewart,  and  A.  Tribble,  2009:  Advances  and  challenges  at  the 
National  Hurricane  Center.  Wea.  Forecasting,  24,  395-419. 

Rotunno,  R.,  and  K.  A.  Emanuel,  1987:  An  air-sea  interaction  theory  for  tropical 

cyclones.  Part  II:  Evolutionary  study  using  a  nonhydro  static  axisymmmetric 
model.  J.  Atmos.  Sci.  44,  542-561. 

Rosenthal,  S.  L.,  1971:  The  response  of  a  tropical  cyclone  model  to  variations  in 

boundary  layer  parameters,  initial  conditions,  lateral  boundary  conditions,  and 
domain  size.  Mon.  Wea.  Rev.,  99,  767-777 . 

Shapiro,  L.  J.,  and  J.  L.  Franklin,  1995:  Potential  vorticity  in  Hurricane  Gloria.  Mon. 

Wea.  Rev.,  123,  1465-1475. 

Smith,  S.  D.,  R.  J.  Anderson,  W.  A.  Oost,  C.  Kraan,  N.  Maat,  J.  DeCosmo,  K.  B. 

Katsaros,  K.  Bumke,  L.  Hasse,  and  H.  M.  Chadwick,  1992:  Sea  surface  wind 
stress  and  drag  coefficients:  The  HEXOS  results.  Bound. -Layer  Meteor.,  60,  109- 
142. 

Smith,  R.  K.,  2006:  Accurate  determination  of  a  balanced  axisymmetric  vortex  in  a 
compressible  atmosphere.  Tellus,  58A,  98-103. 

- ,  M.  T.  Montgomery,  and  S.  Vogl,  2008:A  critique  of  Emanuel’s  hurricane  model 

and  potential  intensity  theory.  Quart.  J.  Roy.  Meteor.  Soc.,  134,  551-561. 

- ,  and  M.  T.  Montgomery,  2009:  Hurricane  boundary  layer  theory.  Quart.  J.  Roy. 

Meteor.  Soc.,  in  review. 

Vickery,  P.  J.,  D.  Wadhera,  M.  D.  Powell,  and  Y.  Chen,  2009:  A  hurricane  boundary 
layer  and  wind  field  model  for  use  in  engineering  applications.  J.  Appl.  Meteor. 
Climat.,  48,  381-405. 

Wang,  J.,  J.  Bian,  W.  O.  Brown,  H.  Cole,  V.  Grubisic,  and  K.  Young,  2009:  Vertical  air 
motion  from  T-REX  radiosonde  and  dropsonde  data.  J.  Atmos,  and  Oceanic 
Technol.,  26,  928-942. 

Zhang,  J.  A.,  P.  G.  Black,  J.  R.  French,  and  W.  M.  Drennan,  2008:  First  direct 

measurements  of  enthalpy  flux  in  the  hurricane  boundary  layer:  The  CBLAST 
results.  Geophys.  Res.  Lett.,  35,  L14813.  doi:10.1029/2008GL034374. 


130 


- ,  W.  M.  Drennan,  P.  G.  Black,  and  J.  R.  French,  2009:  Turbulence  structure  of  the 

hurricane  boundary  layer  between  the  outer  rainbands.  J.  Atmos.  Sci.,  66,  2455- 
2467. 

Zipser,  E.  J.,  R.  J.  Meitin,  and  M.  A.  LeMone,  1981:  Mesoscale  motion  fields  associated 
with  a  slowly  moving  GATE  convective  band.  J.  Atmos.  Sci.,  38,  1725-1750. 


131 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


132 


INITIAL  DISTRIBUTION  LIST 


1.  Defense  Technical  Information  Center 
Ft.  Belvoir,  Virginia 

2.  Dudley  Knox  Library 
Naval  Postgraduate  School 
Monterey,  California 

3.  Michael  Bell 

Naval  Postgraduate  School 
Monterey,  California 

4.  Michael  Montgomery 
Naval  Postgraduate  School 
Monterey,  California 

5.  Kerry  Emanuel 

Massachusetts  Institute  of  Technology 
Cambridge,  Massachusetts 

6.  Patrick  Harr 

Naval  Postgraduate  School 
Monterey,  California 

7.  Russell  Elsberry 

Naval  Postgraduate  School 
Monterey,  California 

8.  Qing  Wang 

Naval  Postgraduate  School 
Monterey,  California 

9.  Francis  Giraldo 

Naval  Postgraduate  School 
Monterey,  California 

10.  Philip  Durkee 

Naval  Postgraduate  School 
Monterey,  California 


133 


